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We review the problem of electron-electron interactions in graphene. Starting from the screening 
of long range interactions in these systems, we discuss the existence of an emerging Dirac liquid of 
Lorentz invariant quasi-particles in the weak coupling regime, and strongly correlated electronic 
states in the strong coupling regime. We also analyze the analogy and connections between the 
many-body problem and the Coulomb impurity problem. The problem of the magnetic instability 
and Kondo effect of impurities and/or adatoms in graphene is also discussed in analogy with 
classical models of many-body effects in ordinary metals. We show that Lorentz invariance plays 
a fundamental role and leads to effects that span the whole spectrum, from the ultraviolet to the 
infrared. The effect of an emerging Lorentz invariance is also discussed in the context of finite 
size and edge effects as well as mesoscopic physics. We also briefly discuss the effects of strong 
magnetic fields in single layers and review some of the main aspects of the many-body problem in 
graphene bilayors. In addition to reviewing the fully understood aspects of the many-body problem 
in graphene, we show that a plethora of interesting issues remain open, both theoretically and 
experimentally, and that the field of graphene research is still exciting and vibrant. 
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I. INTRODUCTION 

One of the most important problems in theoretical 
physics is the understanding of the properties of quantum 
systems with an infinitely large number of interacting de- 
grees of freedom, the so-called many-body problem. In- 
teractions are present in almost all areas of physics: soft 
and hard condensed matter, field theory, atomic physics, 
quantum chemistry, nuclear physics, astrophysics, and so 
on. Interactions between particles are responsible for a 
plethora of effects and many-body states, from the band 
structure of crystals to superconductivity in metals, from 
the quark-gluon plasma in heavy ion collisions to asymp- 
totic freedom in quantum chromodynamics (QCD). It is 
the competition between the kinetic energy of the parti- 
cles, that is, their inertia, and interactions among them 
that leads to the richness and complexity of these differ- 
ent phases. For these reasons, many-body interactions 
are very specific, and the hardest to describe theoreti- 
cally. 

One of the greatest theoretical achievements of the 
last century, the Landau theory of the Fermi liquid 
(Baym and Pethick, 1991), asserts something very sim- 
ple but, at the same time, very deep: that the excita- 
tions of a large (indeed, infinite) collection of strongly 
interacting particles can be described as an equally large 
collection of weakly-interacting quasi-particles that carry 
the same quantum numbers as the original particles. This 
statement is far from trivial. Consider, for instance, the 
behavior of electrons in a metal. The electrons inter- 
act among themselves and with the ions in the crystal 
via strong long-range Coulomb interactions. It is not 
at all clear what is the outcome of this complex inter- 
acting problem. Without having any deep theoretical re- 
sources to treat this problem, except an extraordinary in- 
tuition, visionaries like Paul Drude (Drude, 1900a,b) and 
Arnold Sommerfeld (Hoddeson et ai, 1987) settled the 
foundations for the understanding of this complex prob- 
lem by postulating, shamelessly, that (1) electrons prop- 
agate freely in a non-relativistic (Galilean invariant) way 
(Drude's contribution), and (2) electrons obey Fcrmi- 
Dirac statistics (Sommcrfeld's contribution). Galilean 
invariance dictates that the electrons have a kinetic en- 



where p is the electron momentum and m* is a free pa- 
rameter of the theory that, for lack of a better name, 
is called effective mass. Fermi-Dirac statistics implies 
that electrons carry spin 1/2 and that, in the ground 
state, all states with energy below the so-called Fermi 
energy, Ep, are occupied, and all the states above it are 
empty. With these two basic assumptions and simple 
considerations about electron scattering by defects, the 
Drude-Sommerfeld model was capable of describing ex- 
perimental data of several generations of scientists. 

The understanding of why these two assumptions are 
valid for a strongly interacting problem, such as electrons 
in a metal, had to wait for the development of two major 
concepts: (z) the band structure theory that explains that 
the interaction of the electrons with a periodic lattice of 
ions produces states that, as the plane waves described by 
(1.1), are extended over the entire lattice (Bloch, 1928); 
and (m) the theory of screening^ that is, that metals 
are dynamically polarizable materials and that electrons 
act collectively to screen electric fields in their interior 
(Lindhard, 1954). Hence, long range Coulomb interac- 
tions become effectively short ranged and weak enough 
to give substance to Drude's assumptions. In this case 
the effective mass m* reflects the change in the inertia of 
the electron as it moves around in an effective medium. 
Nevertheless, there are situations when these assump- 
tions fail even in crystalline systems, and that is when 
interesting things happen, namely, the free electron pic- 
ture breaks down. 

In fact, there are many instances where the Fermi liq- 
uid ground state becomes unstable. Electrons not only 
interact with static ions but also with their vibrations, 
the phonons. Electron-phonon interactions, in the pres- 
ence of strong screening, can lead to an effective at- 
tractive interaction between electrons producing a catas- 
trophic Fermi surface instability towards a superconduct- 
ing ground state (Tinkham, 1996). Fermi surface insta- 
bilities also happen in special situations in the presence 
of Fermi surface nesting which can lead to charge and 
spin density wave ground states (Gruner, 1994). Crys- 
tals with inner shell electrons, such as transition metals, 
can also have many-body instabilities due to the strong 
local interactions between the electrons, leading to insu- 
lating states with magnetic properties as in the case of 
Mott insulators (Mott, 1949). Another important case of 
Fermi liquid breakdown is when the electron density is 
very low and the screening disappears. 

Notice that in quantum mechanics the momentum of 
the particle relates to its wavelength. A, by p = h/X 
and hence the kinetic energy (1.1) behaves as K = 
Y? I {2m* X^). If the average distance between electrons is 
£ we sec that the average kinetic energy per electron has 
to be of the order Ek ~ T? ri'J'^ (^m*') where = l/f** is 
the average electron density in d spatial dimensions. On 
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the other hand, the Coulomb interaction is given by: 

y(r) = — , (1.2) 

where e is the electron charge, and eg the dielectric con- 
stant of the medium. Notice that the Coulomb energy per 
electron is of the order Ec ~ e^nj'^ jt^. Thus, the ratio 
of Coulomb to kinetic energy is given by = Eq / Ex oc 
(jiQ/ndY^''' , where uq = [m* / (h^ eojY depends only on 
material properties. Therefore, at high electron densities, 
rid'^ riQ, the kinetic energy dominates over the Coulomb 
energy, which can be disregarded, and the Fermi liq- 
uid description is safe. At low densities, Ud <C uq the 
Coulomb energy is dominant and new electronic phases, 
such as ferromagnetism and Wigner crystallization, can 
become stable (Ceperley, 1978). Therefore, the rela- 
tive strength of the kinetic to Coulomb interactions in 
Galilean invariant systems is completely controlled by 
the electron density. Notice that in all the cases dis- 
cussed above the Galilean invariance was kept intact and 
the driving force for the many-body instabilities was the 
enhancement of the Coulomb relative to the kinetic en- 
ergy. 

With the advent of graphcne (Novoselov et ai, 2004a), 
a two dimensional crystal of pure carbon, this picture has 
changed and a new example of Fermi liquid breakdown 
has emerged in a big way. In graphene, due to its peculiar 
lattice structure, the electrons at the Fermi energy are 
described in terms of an effective Lorentz invariant theory 
where the kinetic energy is given by the Dirac dispersion 
(Castro Neto et ai, 2009a) 

Kg^±vf\p\, (1.3) 

where vp is the Fermi-Dirac velocity, and the ± signs 
refer to two linearly dispersing bands. If we take (1.3) at 
face value and reconsider the argument given above on 
the relevance of the Coulomb interactions we reach very 
different conclusions. For one, the form of the Coulomb 
interaction remains the same as in (1.2), since vp is a ma- 
terial's property and hence much smaller than the speed 
of light, c. This means that the photons which mediate 
the Coulomb interaction arc still much faster than the 
electrons and, thus, the electron-electron interaction can 
be considered as instantaneous. Therefore, the Coulomb 
interaction (1.2) actually breaks the Lorentz invariance 
of (1.3). Secondly, because of the linear scaling of the ki- 
netic energy with momentum, we see that the average ki- 
netic energy per electron has to scale like Eq sa hvpm}/'^ 
and consequently the ratio of Coulomb to kinetic energy 
is given by 



Eg eohvp ' 

and is independent of the electronic density n, depend- 
ing only on material properties and environmental con- 
ditions, such as eo. Here, and from now on, we refer to 



graphene's electron density as n. As the electronic prop- 
erties of graphene are sensitive to environmental condi- 
tions, they will be modified by the presence of other lay- 
ers. In fact, as we are going to show, bilayer graphene 
has properties which are rather different than its mono- 
layer counterpart. Furthermore, due to the same peculiar 
dispersion relation, the electronic density of states, p{E), 
vanishes at the Dirac point, p{E) cx and hence 

graphene is a hybrid between an insulator and a metal: 
neutral graphene is not a metal because it has vanishing 
density of states at the Fermi energy, and it is not an 
insulator because it does not have a gap in the spectrum. 
This means that pristine (or lightly doped) graphene can- 
not screen the long range Coulomb interaction in the 
usual (metallic) way, although it is possible to produce 
electronic excitations at vanishingly small energy. This 
state of affairs makes of graphene a unique system from 
the point of view of electron-electron interactions. The 
long-range interactions lead to non-trivial renormaliza- 
tion of the Dirac quasiparticle characteristics near the 
charge neutrality point, and the resulting electronic state 
can be called Dirac liquid, to be distinguished from the 
Fermi liquid behavior at finite chemical potential (away 
from the Dirac point, where conventional screening takes 
place.) 

The unusual relation between kinetic and Coulomb en- 
ergies not only affects the electron-electron interactions, 
but also the interactions of the electrons with charged 
impurities, the so-called Coulomb impurity problem. In 
a metal described by a Galilean invariant theory of the 
form (1.1), screening also makes the interaction with the 
impurity short ranged, and hence the scattering problem 
effectively reduces to the one of a short range impurity. 
In graphene, because of the lack of screening the situ- 
ation is rather different, and one has to face the prob- 
lem of the effect of the long range part of the poten- 
tial. Scattering by long range interactions has a long 
history in physics and it leads to the issue of logarith- 
mic phase shifts (Baym. 1969). In graphene, because 
of its emergent Lorentz invariance, this issue is magni- 
fied. Since Coulomb interactions between electrons and 
electron scattering by Coulomb impurities are closely re- 
lated issues, one expects that many of the anomalies of 
one problem are also reflected in the other. 

Another interesting consequence of the scaling of the 
kinetic energy with momentum is related to the issue of 
electron confinement. If electrons are confined to a re- 
gion of size L the energy of the states is quantized, no 
matter whether the electrons obey Galilean or Lorentz 
invariance. However, the quantization of energy is rather 
different in these two cases. In a Galilean invariant sys- 
tem, like the one described by (1.1) the energy levels 
are spaced as Ai?o oc 1/L^ while in graphene Lorentz 
invariance, (1.3), implies Ai?G oc \/L. Hence, the size 
dependence of the energy levels in sufficiently small sam- 
ples of graphene is rather different than one would find 
in normal metals. Moreover, since the Coulomb energy 
scales like 1 /L we expect Coulomb effects to be stronger 
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in nanoscopic and mesoscopic graphene samples. 



Furthermore, the fact that graphene is a two dimen- 
sional (2D) system has strong consequences for elec- 
tronic motion in the presence of perpendicular magnetic 
fields. Since a perpendicular magnetic field B leads to 
a quantization of the energy in terms of Landau lev- 
els, and the electrons cannot propagate along the di- 
rection of the field, its effect is singular, in the sense 
that the problem has a massive degeneracy. So, strong 
magnetic fields can completely quench the kinetic energy 
of the electrons that become dispersionless. The elec- 
tronic orbits are localized in a region of the size of the 
magnetic length: £_b = ^ he/ (eB). For a Galilean in- 
variant system, such as the one described by (1.1), for 
p « h/is the kinetic energy per electron is of order 
K w hujc B where tuc — V(™*£b) is the cyclotron 
frequency. On the other hand, for graphene, using (1.3), 
one has Eq « hujQ oc \/B where wa = V^vf/^b, which 
is a consequence of the Lorentz invariance. Notice that 
in both cases the Coulomb energy per electron scales like 
Ec oc e^/(eo£_B) (X ^/B. Hence, in a Galilean invariant 
system the Coulomb energy is smaller than the kinetic 
energy at high fields while for Lorentz invariant systems 
they are always comparable. Thus, one expects Coulomb 
interactions to be hugely enhanced in the presence of 
these magnetic fields. In the 2D electron gas (2DEG) 
this unusual state of affairs is what leads to the fractional 
quantum HaU effect (FQHE) (Laughhn, 1983). 



Given all these unusual circumstances, many questions 
come to mind: How does screening of the long range 
Coulomb interaction work in graphene? Can graphene 
be described in terms a Lorentz invariant theory of quasi- 
particles? Is the Coulomb impurity problem in graphene 
the same as in a normal metal? In what circumstances 
is graphene unstable towards many-body ground states? 
Are there quantum phase transitions (Sachdcv, 1999) in 
the phase diagram of graphene? Do magnetic moments 
form in graphene in the same way as they do in normal 
metals? What is the ground state of graphene in high 
magnetic fields? 



The objective of this review is not to cover the basic 
aspects of graphene physics, since this was already cov- 
ered in a recent review (Castro Neto et ai, 2009a), but 
to try to address some of these questions while keeping 
others open. The field of many-body physics will always 
be an open field because a seemingly simple question al- 
ways leads to another question even more profound and 
harder to answer in a definitive way. In many ways, what 
we have done here is to only scratch the surface of this 
rich and important field, and leave open a large number 
of interesting and unexplored problems. 




c) d) 




FIG. 1 (Color online) a) Honeycomb lattice with the two 
sublattices in graphene. The red arrows are nearest neighbor 
vectors, b) Tight-binding spectrum for the tt — tv* bands. The 
horizontal line intersecting the K point corresponds to the 
Fermi level at half-filling, c) Brillouin zone centered around 
the F point, d) Dirac cone resulting from the linearization of 
the tight-binding spectrum around the K points (blue circles). 

II. CHARGE POLARIZATION AND LINEAR SCREENING 

A. Tight-binding spectrum 

In isolated form, carbon has six electrons in the orbital 
configuration ls^2s^2p^. When arranged in the honey- 
comb crystal shown in Fig. 1(a), two electrons remain in 
the core Is orbital, while the other orbitals hybridize, 
forming three sp^ bonds and one Pz orbital. The sp^ 
orbitals form the a band, which contains three localized 
electrons. The bonding configuration among the Pz or- 
bitals of different lattice sites generates a valence band, or 
TT-band; containing one electron, whereas the antibonding 
configuration generates the conduction band (tt*), which 
is empty. 

From a kinetic energy point of view, the electronic sin- 
gle particle dispersion in graphene is essentially defined 
by the hopping of the electrons between nearest neighbor 
carbon sites in the honeycomb lattice. Unlike square or 
triangular lattices, the honeycomb lattice is spanned by 
two different sets of Bravais lattice generators, forming 
a two component basis with one set for each triangular 
sublattice. Defining a label for electrons sitting in each 
of the two sublattices, say A and B, the free hopping 
Hamiltonian of graphene is 

Ho^-tY, HC^zKCRj)] -f h.c. -^^n,(R,), 

(2.1) 

where acr(Ri), b„(Ili) are fermionic operators for sublat- 
tices A and B respectively, n^CRi) is the number op- 
erator, a I labels the spin and (ij) means summa- 
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tion over nearest neighbors. The two energy scales in 
the Hamihonian are t ~ 2.8 eV, which is the hopping 
energy between nearest carbons, and fj,, the chemical 
potential away from half- filling [see Fig. 1(b)]. In a ho- 
mogeneous system, deviations from half-filling (p = 0) 
are routinely induced cither by charge transfer from a 
substrate (Giovannctti et ai, 2008), by application of 
a back gate voltage (Novosclov et ai, 2005, 2004a. b), 
or else by chemical doping (Calandra and Mauri, 
2007; Griineis et ai, 2009; McChesney et ai, 2010; 
Uchoa et ai, 2008b). 

In momentum space the free Hamiltonian of graphene 

is 

p.cr \ P / 

where ^'p.^ = (flp.cr, &p,cr) is a two component spinor and 

3 

0p = ^e'P- (2.3) 

i=l 

is a tight-binding function summed over the nearest 
neighbor vectors 

a \/3 a \/3 

ai^ax, a.2 = --x + a—y, as = --x - a— y , 

(2.4) 

where a « 1.42 A is the carbon-carbon spacing. The di- 
agonalization of Hamiltonian (2.2) yields the spectrum 
of the two TT-bands of graphene in tight-binding approx- 
imation (Wallace, 1947), 

i?±(p) = ±#p|-M- (2.5) 

The +(— ) sign in the spectrum corresponds to the con- 
duction (valence) band. 

The hexagonal Brillouin zone (BZ) of graphene shown 
in Fig. 1(c) has three high symmetry points: the F point, 
located at the center of the BZ, the M point, which 
indicates the position of the Van Hove singularities of 
the TT-TT* bands, where the density of states (DOS) is 
logarithmically divergent, and the K points, where the 
TT-bands touch, and the DOS vanishes linearly. An ex- 
tensive description of the band structure of graphene 
and its electronic properties is reviewed in detail by 
Castro Neto et ai, 2009a. 

B. Dirac fermion Hamiltonian 

The topology of the Fermi surface in undoped 
graphene is defined by the six K points where the con- 
duction and valence bands touch, E±{K.) = ±|0k| = 0. 
These special points form two sets of nonequivalent 
points, K and K' , with K = -K' and |K| = 47r/(3v^a), 
which cannot be connected by the generators of the 
reciprocal lattice. The linearization of the spectrum 



around the valleys centered at ±K gives rise to an 
effective low energy description of the electrons that 
mimics the spectrum of masslcss Dirac particles. In 
this effective theory, the elementary excitations around 
the Fermi surface are described by a Dirac Hamiltonian 
(Semenoff, 1984), 

•Ho = 51 *L • 7 - htq (g) (Jo] *k<T, (2.6) 

ak 

where 

^ko- = (aK+k,(T, &K+k,(T, fe-K+k,o-, fl-K+k.o-) (2.7) 

is a four component spinor for sublattice and valley de- 
grees of freedom. In this representation, 7^ = T3 (g) ai 
, where r and <t are the usual Pauli matrices, which 
operate in the valley and sublattice spaces respectively 
(i = 1,2,3 correspond to x, y and z directions, and tq = 1 
and ctq = 1 are identity matrices) . The form of the spec- 
trum mimics the relativistic cone for massless fermions 
(Wallace, 1947), 

E±{k) ^ ±v\k\ - fi (2.8) 

where the Fermi velocity v — (3/2)ta « 6eVA is nearly 
300 times smaller than the speed of light, i.e. u w 1 x 
lO^m/s. From now on we set h ^ ks = I everywhere, 
except where it is needed. For simplicity of notation, we 
call the Fermi velocity v (i.e. vp = v) throughout this 
review. 

The Hamiltonian (2.6) is invariant under a pseudo- 
time reversal symmetry operation, S = i(to <8) cr2)C, 
SHS~^ = H, (C is the complex conjugation operator), 
which is equivalent to a time reversal operation for each 
valley separately. It is also invariant under a true time 
reversal symmetry (TRS) operation, which involves an 
additional exchange between the valleys, T = (n (g) <7i)C. 

In the absence of back scattering connecting the two 
valleys, the Hamiltonian can be decomposed in two in- 
dependent valley species of Dirac fermions with opposite 
chiralities: 

Ho,+ = J2^l ^^Jvk-cT- (2.9) 

(r,k 

Ho,- = ;^*L^kJ-^'k-'T* -A*]*-,k., (2.10) 

o-,k 

where ^'±.kcr = (a±K+k.(T, ^±K+k.cr) arc two component 
spinors. In this review, unless otherwise specified, we 
will arbitrarily choose one of the two cones and assume 
an additional valley degeneracy in the Hamiltonian. So 
valley indexes will be generically omitted unless explicitly 
mentioned. A more detailed description of the symmetry 
properties of the graphene Hamiltonian can be found in 
(Gusynin et ai, 2007). 



FIG. 2 Diagram for the polarization bubble corresponding 
to eq. (2.12). 



C. Polarization function 

The Green's function of graphene is a 2 x 2 matrix 
represented in the sublattice basis by 



G(k,r) 



Gaa Gab 
Gba Gbb 



where Gaa = ~ {T[a'k.{T)a\^{0)\) and so on, with r as the 
imaginary time. In the low energy sector of the spectrum, 
close to the Dirac points, the non-interacting Green's 
function is G^*^^ (k, iw) = [zw + /i — uk • cr] ^. or equiv- 
alently, in a chiral representation. 



1 + S(Tk 



s=± 



fj, — sv\\i\ 



(2.11) 



where = cr • k/|k| is twice the quantum mechanical 
helicity operator for a Dirac fermion with momentum k, 
and s = ± labels the two branches with positive and 
negative energy in one cone. It is clear that the positive 
and negative branches within the same cone have also 
opposite helicities. 

The polarization function in one loop is calculated di- 
rectly from the bubble diagram shown in Fig. 2, 



n(i)(q,H-^EE-^-'.-'(P'q)^ 

f[Es,{p + <l)]- f[Es{p)] 
Es'{p + q) - Es{p) ~ iuj 



(2.12) 



where f{E) = (e^/-^ + l) is the Dirac-Fermi distribu- 
tion, with T as temperature, iV = 4 is the degeneracy for 
two spins and two valleys, and 



J"s,s/(p,q) = itr(l 



(2.13) 



are the matrix elements due to the overlap of wavefunc- 
tions for intraband (s = s') and interband (s = — s') tran- 
sitions, 'tr' means trace over the sublattice indexes. In 
a more exphcit form, J^s,s'{p, q) = [1 + ss' cos6'p_p+q] /2, 
where 9 is the angle between p and p + q. The full mo- 
mentum, frequency, and chemical potential dependence 
of (2.12) is shown in panels (a-d) of Fig. 3. 

In metals, screening is a many-body property directly 
related to the polarizability of the electrons around the 
Fermi surface. In graphene, because the density of 
states (DOS) vanishes linearly around the Dirac points, 
p{E) (X \E—ii\/v'^, exactly at the neutrality point (/i = 0) 



(a) Ren{q,cj)/p(^i) (b) Im n(q, w) / p(p) 

■3 4f| 




2 3 4 5 

q/kp 





Im 


Re / 



FIG. 3 (Color online) Polarization bubble Il'^-'(g,cj) for 
graphene, within the Dirac approximation. Panels (a) and 
(b) show, respectively, a density plot of the real and imagi- 
nary parts of the polarization bubble, Tl''-^\q,uj), defined in 
eq. (2.12), and normalized to the DOS at the Fermi level, 
p(/i). Panels (c) and (d) present constant frequency cuts at 
dj/^i = 0.5, 1.0, 1.5, 2.0, 2.5, 3.0. In panel (e) we show the 
static limit, n''^^(q', 0), whose closed form expression is writ- 
ten in eq. (2.16). Notice the transition from a constant value 
{q < 2kF) to the linear in q dependence at large momenta. 
The derivative of the polarization is shown in the same panel, 
and can be seen to vary continuously. In (f) we plot the real 
(black/solid) and imaginary (red/dashed) parts of the uni- 
form limit (2.17). 



the screening of charge is completely suppressed, and the 
polarization function describes the susceptibility of the 
vacuum to particle-hole pair production, exactly as in 
the diagonal time component of the polarization tensor 
in massless Quantum Electrodynamics (QED), QED2+1 
(Appelquist et ai, 1988; Gonzalez et ai, 1994; Pisarski, 
1984), 



n(i)(g,c 



1 



(2.14) 



4 y/ v^q'^ — 

Here we have performed a Wick rotation to real frequen- 
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cies, iw — ^ w + 0+. Since the Fermi surface in this case 
is just a point, there is no phase space for intraband ex- 
citations at zero temperature due to the Pauh principle. 
The process of creation of particle-hole pairs involves in- 
coherent excitations of electrons from the lower to the 
upper band. The continuum of particle-hole excitations 
is well defined for all virtual transitions with w > vq. 

For finite fi there is a crossover in the behavior of the 
polarization function. The DOS around the Fermi level is 
finite and the intraband excitations dominate the infrared 
behavior of the polarization. For vq <^ \^\ and w <C \fj,\, 
the leading term in the polarization function is (Shung, 
1986a) 

As in a Fermi liquid, there is a particle-hole continuum 



for uj < vq, which is due only to intraband transitions. 
The polarization function in graphene is a regular func- 
tion everywhere except at \uj\ ~ vq, where it has an on- 
shell singularity delimiting the border of the particle-hole 
continuum. 

The polarization was derived originally by Shung, 
1986a and later rederived by a number of authors (Ando, 
2006; Barlas et ai, 2007; Hwang and Das Sarma, 2007; 
Wunsch et ai, 2007). These results rely on the cone ap- 
proximation, which ignores contributions coming from 
the non linear part of the spectrum. In addition, the band 
width is assumed to be infinite. Although the charge po- 
larization for Dirac fermions in 2D is well behaved and 
does not require cut-off regularization in the ultraviolet, 
the physical cut-off of the band, D, generates small cor- 
rections that vanish only in the _D — ^ od limit. In this 
sense, the 'exact' expression for the static polarization 
function (uj = 0) for arbitrary momentum is 



nW(g,0) = -^ + %-2M;^ 
nv Zttv 



2kF 



•l- 



'2kF'' 



TT 

2 



(2.16) 



where kp = \^i\/v is the Fermi momentum, and Oix) 
is a step function. The static polarization is plotted in 
Fig. 3(e). 

At qK, 2kF the static polarization exhibits a crossover 
from a two dimensional electron gas (2DEG) to Dirac 
fermion behavior. For details of the polarization func- 
tion in the 2DEG please refer to Fig. 4. As in the 2DEG, 
the polarization of graphene is constant for q < 2kp. 
For q > 2k F, it eventually becomes linear in q for large 
momenta. At the crossover, the static polarization and 
its first derivative are continuous at q — 2k p. The dis- 
continuity only appears in the second derivative. This is 
distinct from the 2DEG case, where the first derivative is 
discontinuous. The difference will affect the spacial de- 
pendence of the Friedel oscillations in the two systems. 

In the opposite limit, for arbitrary uj and q — >■ 0, the 
polarization function becomes 

n«(<z-o,.) = -^[^ + iinffc^y , 

(2.17) 

which is shown in Fig. 3(f). The presence of a pocket 
of electrons (holes) around the Dirac points opens a gap 
in the particle hole continuum for interband excitations 
(w > vq). From Eq. (2.17), it is clear that the imagi- 
nary part of the polarization function at small momen- 
tum is zero unless uj > 2\^\ [Fig. 3(b)]. This is so be- 
cause the phase space for vertical interband excitations 
is Pauli blocked for w < 2|/i|, generating a gap for optical 
absorption in the infrared. At finite q, the threshold for 



interband transitions is w > 2\fi\ — vq for q < 2kp, as 
shown schematically in Fig. 5. 



D. Collective modes and screening 

The Coulomb interaction among the electrons in 
graphene gives rise to collective modes and metallic 
screening when the Fermi level is shifted away from the 
Dirac points. In a 2D system, the bare Coulomb interac- 
tion is given by 

V{q)^^, (2.18) 

where e is the charge of the electron and eo is the effec- 
tive dielectric constant of the medium. For graphene in 
contact with air and a substrate with dielectric constant 
K, eo = (1 -l- k)/2. In most of the experiments, graphene 
lies on top of some substrate like Si02 or SiC, where di- 
electric effects are moderate (for instance, the dielectric 
constant of Si02 is k w 4). The background dielectric 
constant can be significantly enhanced in the presence of 
substrates in contact with strong dielectric liquids such 
as ethanol (k « 25) or water (k sa 80) (Jang et ai, 2008; 
Ponomarenko et ai, 2009). 

As usual, the collective modes follow from the zeros of 
the dielectric function 

eiq,Uj)=eo[l~Viq)n^'Hq,uj)], (2.19) 
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Re n(q, oj)lp(n) 
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FIG. 4 (Color online) Polarization bubble n''^'(g,aj) for the 
conventional 2DEG. Panels (a) and (b) show, respectively, a 
density plot of the real and imaginary parts of the polarization 
bubble, n'^' (5, a;), normalized to the DOS at the Fermi level. 
Panels (c) and (d) present constant frequency cuts at uj/ ji = 
0.5, 1.0, 1.5, 2.0, 2.5, 3.0. In panel (e) we plot the static limit, 
n^^'(g,0), and in (f) the uniform limit, n(^)(0,aj). 



calculated here in the Random Phase Approximation 
(RPA). Since graphene is a 2D system, the collective 
plasmon mode is gapless. The leading term in the po- 
larization for small frequency and momenta (compared 
to kp) is shown in Eq. (2.15). From it one can easily 
extract the infrared dependence of the plasmon. 



ujp{q) = v^(2^eVeo)g , 



(2.20) 



which follows the same dispersion as the plasmon en- 
countered in the 2DEG. The ^ dependence of the plas- 
mon was recently confirmed by a high resolution en- 
ergy loss spectroscopy (EELS) measurement in graphene 
(Liu et ai, 2008). Additional corrections due to the in- 
terband excitations (which are absent in the 2DEG) can 
be absorbed into the definition of the background diclcc- 




00 
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FIG. 5 (Color online) Colored regions represent the particle- 
hole continuum of graphene due to interband (gray area) and 
intraband (green) transitions. On the left: half-filled case; 
right: finite jj, case, away from half filling. Dashed line: acous- 
tic plasmon for the single layer (cjp oc ^/Jlq). 



trie constant (Shung, 1986a), 



2u}p{q) \2\^\+u}p{q) 



(2.21) 



As in the 2DEG, the screened Coulomb interaction for 
q < 2kp is 



V{q) 



1 



e((7,0) eoq + qTF 



(2.22) 



where qrp = ^Tre^kp / [veo) is the Thomas-Fermi momen- 
tum [kp = \ii\lv), which sets the size of the screening 
cloud. In the presence of an external charged impurity 
Ze, the induced charge, 5Z, has a non-oscillatory com- 
ponent coming from the g — )■ limit of the polarization 
that decays as {kpr^)~^ (as in a 2DEG), and an oscilla- 
tory part which corresponds to the Friedel oscillations at 
q = 2k p. The Friedel oscillations in graphene decay as 
cos{2kpr) / (kpr^) , differently from the 2DEG case, where 
the decay is of the form cos{2kpr)/r'^ . The difference is 
caused by the fact that the static polarization function in 
the 2DEG has a cusp at g = 2kp, whereas in graphene, 
the first derivative is continuous [cfr. Figs. 3(e) and 4(e)]. 

For undoped graphene, V{q)Il^^'> = — (7r/2)[e^/(weo)] 
[see Eq. (2.14)], and the static dielectric function is a 
constant. The effective Coulomb interaction in this case 



IS 



e(<Z,0) 



1 27re^ 
eRPA q 



(2.23) 



where erpa = eo -f (7r/2)(e^/u) is the effective back- 
ground dielectric constant, renormalized by the inter- 
band transitions. Additional many body effects resulting 
from self-energy insertions in the bubbles logarithmically 
renormalize this correction to zero in the g — >■ limit, 
as will be clear in Sec. Ill of this review. On the dy- 
namical side, inserting Eq. (2.14) into Eq. (2.19), one 
can easily see that no collective modes are allowed in un- 
doped graphene, at zero temperature, within the RPA 
framework. At half-filling, RPA is justified in the limit 
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of large number of fermionic species, N, which favors di- 
agrams with maximal number of bubbles at each order of 
perturbation theory. In graphene, the physical number 
of species is iV = 4, and additional corrections beyond 
RPA coming from the exciton channel near the on-shell 
singularity of the bubble, |a;| ~ vq, were shown to gener- 
ate a new acoustic plasmon mode (Gangadharaiah et al., 
2008). In the static limit (w — )■ 0), vertex corrections in 
the bubble are perturbatively small and RPA can be jus- 
tified in the calculation of the dielectric function even at 
half- filling (Kotov et ai, 2008b). The structure of per- 
turbation theory in graphene will be discussed in detail 
in Sec. III. 

In addition to the low energy acoustic mode due to 
intraband transitions, graphene has also two high en- 
ergy optical plasmons generated by interband excitations 
around the Van-Hove singularities of the tt — n* bands, 
and also by optical transitions between u — tt* and n — a* 
bands (Eberlein et ai, 2008; Kramberger et ai, 2008). 
The measured optical gaps of the tt and tt — a band 
plasmons in graphene are 4.5 eV and 15 eV, respectively. 
Similar modes were also observed in graphite, where they 
appear blue shifted to 7eV and 24 eV respectively, ac- 
cording to optical data (Taft and Philipp, 1965), X-ray 
measurements (Shulke et ai, 1988), and ab-initio calcu- 
lations (Marinopoulus et at, 2004). 



E. Infinite stack of layers 

In the case of an infinite stack of graphene layers, the 
Hamiltonian term for the Coulomb interaction among all 
the electrons can be written in real space as 

p2 r 1 
'Hc = — d^rd^r' n(r)^ -n(r') , (2.24) 

where n(r) is the 3D particle density operator. In the 
absence of interlayer hopping, as in the case for exam- 
ple of several graphite intercalated compounds, the elec- 
trons remain confined in each layer, but the unscreened 
Coulomb lines fill the entire space in between the layers, 
coupling all the electrons in the system. In that case we 
may constrain the local density operator h to be in the 
form (Visscher and Falikov. 1970) 



n(r) 



OO 

/— — oo 



h{r)d{z - Id) 



(2.25) 



where / is an integer labeling the layers, and d is the 
distance between layers. In momentum space, making 
a discrete sum over the layers, the Coulomb interaction 
between all the electrons is 



where 



l/(k) =27rd S{q,k,) 



(2.26) 



(2.27) 



with k = (q, fcz), q is an in-plane momentum, and 
(Fetter, 1974) 



sinh(qd) 



co&]i{qd) — cos(kzd) 



(2.28) 



is the structure factor for a stack with an infinite num- 
ber of layers. In the limit when the distance between the 
layers d is small, Eq. (2.27) recovers the isotropic case 
V{k) = 47r(e^/eo)/(9^ + fcz), whereas in the opposite limit 
[d oo) one gets the 2D case, V{k) = 2tt d{e'^ /eo j/q. 
In any case, the polarization function must be integrated 
over a cylindrical Fermi surface of height iir/d, and so 
, w) acquires an additional factor of l/d compared 
to the single layer case. The extension of this problem to 
include the interlayer hopping dispersion in the polariza- 
tion was considered by Guinea, 2007. 

Away from the neutrality point (/x ^ 0), instead of 
a single acoustic mode as in the monolayer, the ze- 
roes in the dielectric function of the multilayer gener- 
ate a plasmon band, where the modes are labeled by 
kz G [— tt/c?, Tr/d]. For q <C l/d, the plasmon dispersion 
is (Shung, 1986a) 



2/xe= 



qS{q,kz 



(2.29) 



In the fcz = mode, the charge fluctuations be- 
tween different layers are in-phase, and the result- 
ing plasmon mode is optical, ajp(g,0) ~ (A^e'^/ead) + 
|(ug)^. For UJp{q) > 2fi, this mode is damped by 
the particle-hole continuum due to interband transitions 
(see Fig. 6), in agreement with energy loss spectroscopy 
data (Laitenberger and Palmer, 1996). The out-of-phase 
modes (for k^ ^ 0) arc acoustic. At the edge of the plas- 
mon band, the mode k^ = ±7r/d disperses linearly with 
the in-plane momentum, LOp{q^±'K / d) = fie-^d/eoq, in 
contrast with the 2DEG dispersion (wp oc y^) present 
in the single layer. Except for the lack of an inter- 
band particle-hole continuum and the associated damp- 
ing, similar plasmon band features are also expected in 
the 2D layered electron gas, for fermions with quadratic 
dispersion (Hawrylak, 1987). 



F. f-sum rule 

The /-sum rule is a generic statement about conser- 
vation of the number of particles and results from the 
analytical properties of the retarded charge susceptibil- 
ity. It can be generically defined as (Nozieres, 1964) 



dujujlmx^(k,uj) ^TT{[[n,n{-k)],h(k)]), (2.30) 



where Ji is the Hamiltonian, h is the particle density op- 
erator, is a retarded charge susceptibility, x(k, t) ~ 
(r[n(k, T)n(— k, 0)]), and (...) is an expectation value cal- 
culated in some basis. 
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FIG. 6 (Color online) Plasmon band (hatched region) for 
an infinite stack of graphene layers. Red line: optical mode 
kz = 0. Dashed line: acoustic mode kz = Tr/d, Up oc y^q, 
with linear dispersion, at the edge of the band. All the other 
modes in between are acoustic. Adapted from Shung, 1986a. 

As in any solid, the exact electronic Haniiltonian of 
graphene can be decomposed into a Hamiltonian of free 
electrons, plus a periodic potential due to the lattice, and 
interactions. If the interactions depend only on densities, 
the commutators in Eq. (2.30) can be calculated exactly, 
and the only term that survives is the Kinetic energy due 
to the free electrons, 

{[[n,h{~k)],h{k)])^N,— , (2.31) 

m 

where m is the bare electron mass and Ng is the number 
of fermions in the band. Choosing, for example, a basis 
of non-interacting fermions, the sum rule in graphene is 




duj a;Imn(i) (fc, w) = tt^— , (2.32) 
m 



as in metals, where n'^''(fc,a;) is the bare polarization 
bubble, calculated using the full non-interacting spec- 
trum (dictated by the lattice symmetry). The validity 
of the / sum-rule does not require Galilean invariance of 
the quasiparticles, but of the free electrons, which are not 
relativistic and hence obey the Schrodinger equation. 

For low energy effective Hamiltonians, such as the 
Dirac Hamiltonian in graphene (which do not include 
the periodicity of the spectrum in the Brillouin zone), 
the /-sum rule above is still formally satisfied when ap- 
plied for the electrons (holes) in the conduction (valence) 
band only, as can be explicitly checked by direct substi- 
tution of the polarization due to intraband transitions, 
Eq. (2.15), into Eq. (2.32). The number of electrons 
(holes) in this band, Ne = kpA/n, where A = 3%/3a^/2 
is the unit cell area, is set by the size of the Fermi surface, 
and the verification of the sum rule follows as in a Fermi 
liquid. 

The Dirac Hamiltonian, however, violates the /-sum 
rule (2.32) when interband transitions are taken into ac- 
count. In that case, the left hand side of Eq. (2.32) be- 
comes independent of the chemical potential, consistent 
with the fact that (Sabio et ai, 2008) 

m,hi-k)],h{k)])^k'^ (2.33) 



for a Dirac Hamiltonian, where D is the ultraviolet cut- 
off. A similar dependence with the cut-off also occurs 
in the true 3D relativistic problem, where the sum rule 
reflects the number of particles contained in the vacuum 
of the theory, which is formally divergent (Ceni, 2001; 
Goldman and Drake, 1982). In graphene, as in any two 
band semi-metal or semiconductor, the validity of the /- 
sum rule is physically recovered when the periodicity of 
the electronic spectrum is restored back into the Hamil- 
tonian. 



III. QUASIPARTICLES IN GRAPHENE 

The quasiparticle properties of graphene are modi- 
fied by the presence of long-range Coulomb interactions. 
Their effects are especially pronounced when the Fermi 
energy is close to the Dirac point (^ sa 0), and can re- 
sult in strong renormalization of the Dirac band struc- 
ture (the Fermi velocity v), and the quasiparticle residue 
(Z). Consequently, many physical characteristics, such 
as the compressibility, spin susceptibility and the spe- 
cific heat can be strongly affected by interactions. Even 
when the Fermi surface is large and the system is a Fermi 
liquid, there arc strong modifications of the physics near 
the Dirac point due to the presence of additional peaks in 
the quasiparticle decay rate, related to plasmon-mediated 
decay channels. Even reconstruction of the Dirac cone 
structure near the charge neutrality point appears possi- 
ble, as indicated by recent Angle- Resolved Photoemission 
Spectroscopy (ARPES) measurements. All these effects 
are sensitive to the value of the Coulomb interaction con- 
stant in graphene, a. 

A. Low-energy behavior near the Dirac point 

1. Weak-coupling analysis 

The interaction parameter which characterizes the 
strength of the Coulomb interaction in graphene is 
(Eq. (1.4)) 



At kp =0 screening is absent, and the interaction po- 
tential in momentum space: 

27re2 

Vip) = (3.2) 

The value of a ~ 2.2/eo depends on the dielectric envi- 
ronment since, as previously discussed, eo = (1 + k)/2 
for graphene in contact with air and a substrate with 
dielectric constant k. In vacuum, a ~ 2.2. 

In the case of small coupling, a <C 1, we can employ 
standard perturbation theory, involving the perturbative 
computation of the self-energy I](fe, uj), which enters in a 
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FIG. 7 Self-energy diagrams: (a) First order Hartree-Fock, 
(b) Second order loop diagram (first diagram in the RPA se- 
ries), (c) Second order exchange (vertex correction) diagram, 
(d) Rainbow diagram. 



standard way the Dirac fermion Green's function (GF), 
for a given valley: 



G{k,uj) 



1 



WCTo ^ vcr ■ k — I](fc,a;) -|- i(7oO+sign(a;) 



(3.3) 



It is convenient to decompose the self-energy into two 
pieces with different pseudo-spin structure 

I](fc, oj) = I]o(fc, w) + St,(A;, a;), Eq oc croi St, (x <t • k, 

(3.4) 

where ctq = 1 is the unit matrix, which from now on will 
not be written explicitly. Then we have 



G{k,uj) 



uj — Zivcr ■ k + Et,) 
where Z is the quasiparticle residue 

Z"^ = 1 - dJ^o/duj, 



(3.5) 



(3.6) 



and E„ is responsible solely for the velocity renormaliza- 
tion. 

The first order diagram shown in Fig. 7(a) is the 
Hartree-Fock exchange contribution, and can be readily 
evaluated (we denote by G'-^' the non-interacting GF): 

sW(fc,..)=* I |^GW(fc+p,.. + e)np), (3.7) 

which at low external momenta exhibits a logarithmic 
singularity 

I](i)(fc,a;) = Y^iPik) = ^v(t -kMA/k), A/k > 1. 

(3.8) 

At this order we have Eo = 0, i.e. Z = 1 due to the 
frequency independence of the interaction potential, and 
the quasiparticle velocity increases: 

u(fc) = w (^1 + |ln(A/fc)) , A/fc>l. (3.9) 

The ultraviolet cutoff A ~ 1/a represents the momentum 
scale up to which the spectrum is Dirac-like. 



FIG. 8 (a) Self-energy and (b) Vertex corrections to the po- 
larization bubble 



While the linearity of the spectrum in graphene was 
realized a long time ago (Wallace, 1947), in the con- 
text of studying graphite formed by layers of graphene, 
the self-energy correction Eq. (3.8) due to interac- 
tions was first investigated perturbatively much later 
by Gonzalez et ai, 1994. The non-trivial velocity 
renormalization is due to the unscreened, long-range 
Coulomb interactions. Similar logarithmic divergencies 
were also found in gaplcss 3D semiconductors, where 
the Dirac spectrum originated from special symmetries 
(Abrikosov and Beneslavskii, 1971). 

The above calculation forms the basis of the Renor- 
malization Group (RG) analysis. In the RG spirit one 
integrates out the high momentum degrees of freedom, 
i.e. regions of momenta A > |p| > Ai, and the results 
vary with the quantity ln(A/Ai) = I. Here we denote by I 
the RG parameter, so that the infrared limit corresponds 
to / — > oo (i.e. one integrates down to the infrared scale 
k^O,l^ ln(A/fc)). From Eq. (3.9) we obtain 



dv 
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e 

4eo 



(3.10) 



This equation has to be supplemented with an additional 
equation reficcting the absence of charge (e^) renormal- 
ization: 



de2 
'df 







(3.11) 



There are several ways to understand this. It was ar- 
gued early on that the vertex function does not acquire 
any divergent contributions, which is related to the ex- 
pected regular behavior of the polarization operator to 
all orders in graphene (Gonzalez et al, 1994). More re- 
cently, explicit calculations up to two loop order were 
performed (de Juan et ai, 2010; Kotov et ai, 2008b); it 
was confirmed that the vertex function is finite in the low- 
energy limit. In addition, direct examination of the po- 
larization function at two loop level (Kotov et ai, 2008b) 
found that the self-energy correction. Fig. 8(a), acquires 
a logarithmic divergence which can be absorbed into the 
renormalized velocity v{k) (Eq. (3.9)), while the vertex 
correction of Fig. 8(b) is finite: 



27re^ 



H(2f')(q,0) = finite = -O.SSa^ 



(3.12) 



Incidentally, this contribution leads to enhancement of 
the dielectric static screening (i.e. the dielectric constant 
beyond linear (RPA) order becomes e = 1 



fa-hO.53 a^.) 
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Alternatively, one can argue that in two-dimensional 
field theories with Coulomb interactions the charge 
does not flow because it appears as a coefficient 
in a nonanalytic term in the action (Herbut, 2006; 
Ye and Sachdev, 1998). The conclusion then is that only 
the quasiparticlc velocity and residue (see below) are 
renormalized. In particular, at first order we can combine 
Eqs. (3. 10), (3. 11) into a single one refiecting the renor- 
malization (running) of the coupling a: 



da 



a 



(3.13) 



Therefore we have an infrared stable fixed point at a = 0, 
and the flow towards it is logarithmic: 



a{k) 



ln(A/fc)' 







(3.14) 



Thus the Coulomb interactions are marginally irrelevant. 
This is equivalent to a logarithmically divergent velocity: 
v{k) - (eV4)ln(A/fc), fc ^ 0. 



a. Two-loop results. It is instructive to examine 
corrections beyond first order (Mishchenko, 2007; 
Vafek and Case, 2008), since additional effects appear, 
such as renormalization of Z. For example the first 
diagram in the RPA series shown in Fig. 7(b) is 

^^''Hk,c.)^^ J |^G("'(fc+P,-+e)(np))^n(i)(p,£). 

(3.15) 

An explicit evaluation at low energies and momenta gives 
a single logarithmic divergence 



I](2f)(fc,c.) 



"24" 



{uj + v<T- k) ln(A/A:), fc/A 0, 



I.e. 



(2b) _ 



(3.16) 

-(A^aV24)(a;)ln(A/fc), and sl,^''^ = 
— {Na^ /2A)va ■ k ln(A/fc). Because the polarization bub- 
ble is proportional to the number of fermion flavors 
A'^ = 4 (valley -l-spin), we have explicitly written the 
N dependence. By comparing with Eq. (3.5), we find 
that the velocity is changed by an amount {—N/24: — 
N/24)a^v\n (A/fc). 

In addition, other diagrams at second order have to be 
added, such as the vertex correction of Fig. 7(c). Most 
importantly, this diagram is also proportional to In A. 
Collecting all contributions one finds the RG equation 
for the velocity flow (Vafek and Case, 2008) 



dv 



a 



N 
12 



(3.17) 



with S w 0.03. One observes that the contribution of the 
"RPA" diagram is numerically dominant at second order 
(it is larger than the rest by a factor of 10 for N = 4.) 
In addition, the second order tendency is a decrease of 
the velocity. Consequently a flnite coupling fixed point is 



possible at ac « 0.8. This fixed point is infrared unstable 
since near ^ = —C{a — ac)v, C > 0, i.e. for a > 

V flows towards zero {a flows to oo) while for a < ac, 

V flows towards oo (a flows to zero.) Of course it is not 
clear that this estimate is reliable since the fixed point 
value ac is not small, and we used perturbation theory 
(a ^ 1) to derive this result. On the other hand, a flow 
towards strong coupling for a > ac is consistent with 
the formation of an excitonic insulator (mass generation) , 
for which strong evidence has accumulated by now, as 
we discuss in Section III.B. Recent numerical simulations 
give the value ac ~ 1 (see Section III.B). 

Finally, we also find that Z is renormalized at second 
order, since the self-energy is frequency dependent. From 
Eq. (3.6) we can expand to second order of bare pertur- 
bation theory Z « 1 — ln(A/A:), which would lead us 

to an RG equation for Z: ^ — -~^^Z, to be solved 
together with Eq. (3.13), or Eq. (3.17), depending on the 
desired level of approximation. Alternatively, Eq. (3.6) 
is already written in a "nonperturbative" way. Ignoring 
for the moment the running of a, we have at low energies 



Z = 



1 



24 



l + i^ln(A/fc) NaHniA/k)' 



fc/A ^ 0. 



(3.18) 

This result, along with the previous one for Sg^''^ brings 
us to the infrared behavior (we use co and fc interchange- 
ably in the infrared limit): 

So - a2w|ln(w)|, w^O. (3.19) 



a2|ln(w)|' 



This is characteristic of a marginal Fermi liquid 
(Das Sarma et ai, 2007; Gonzalez et al., 1994). How- 
ever, this regime is never achieved if the running of a is 
taken into account, as is intuitively clear from the above 
equations. As we will see later from the solution of the 
RG equations for Z and a, in fact Z tends to level off in 
the infrared, and the system has well-defined quasiparti- 
cles. 

It is interesting to note that trigonal distortions, 
which change the band structure away from the Dirac 
equation, are modified by the electron-electron interac- 
tion, and their irrelevance at low energies is enhanced 
(Foster and Aleiner, 2008). As a result, the linear dis- 
persion becomes an even more robust feature of graphene 
(Roldan et al, 2008). 



b. Influence of disorder. Before we proceed, let us briefly 
address the effect of disorder. Two major sources of 
disorder are scalar potential random fluctuations (e.g. 
formation of electron- hole puddles), and vector gauge 
field randomness, related to formation of ripples. Start- 
ing with the latter, i.e. a gauge field coupled to the 
Dirac fermion pseudospin cr • A, and characterized by 
variance A, (A^(ri)^i,(r2)) = ASf^^d{ri — r2), one can 
readily derive the corresponding RG equations in the 
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Disorder A 



FIG. 9 An attractive line of fixed pints for interactions and 
gauge field disorder. 



weak disorder and interactions limit (Herbut et at, 2008; 
Stauber et ai, 2005) 

dA ^ da A 

Gauge field disorder itself is not renormalized, while the 
interplay of disorder and interactions leads to a line of 
attractive fixed points located at: a* ~ ^A, as shown 
in Fig. 9. Physically the variance is related to the char- 
acteristic height h, and length L of the corrugations of 
the surface, A ~ h'^/{L'^a^). Thus weak disorder generi- 
cally shifts the fixed point away from a = 0, while strong 
disorder can have an even more profound effect (Section 
VI.C). 

In addition, for weak interactions, the inclusion of 
scalar (density fluctuations) disorder turns out to be 
a relevant perturbation which grows under renormal- 
ization, and thus away from the perturbative regime 
(Aleiner and Efetov, 2006). Moreover, gauge field disor- 
der, when combined with strong-enough interactions, can 
cause the interactions to grow (Vafek and Case, 2008). It 
has been argued that the strong-coupling regime for dis- 
order and interactions generically occurs when all types 
of disorder consistent with graphene's symmetries are in- 
cluded (Foster and Aleiner, 2008). 

A detailed analysis of this complex situation is beyond 
the scope of this work, and from now on we continue our 
discussion of clean graphene. 



2. Strong-coupling/RPA analysis 

The full RPA treatment was performed by many au- 
thors (Das Sarma et ai, 2007; Foster and Aleiner, 2008; 
Gonzalez et ai, 1999; Kotov et ai, 2009; Polini et al., 
2007; Son, 2007). Here we mostly follow Son, 2007. The 
RPA self-energy is shown diagrammatically in Fig. 10, 
and corresponds to the equation 
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FIG. 10 RPA self-energy, which includes an infinite resum- 
mation of polarization bubbles. 



The RPA potential is given by 



eop- 27re2n(i)(p,e)' 



(3.22) 



Quite remarkably, at low momenta one can evaluate the 
singular contribution to the self-energy analytically 



where we have defined 



-Fo{\)uj + Fi{\)vcT ■ k] ln(A/fc) , 

(3.23) 
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8 



-Na. 



(3.24) 



This parameter is measuring the importance of polariza- 
tion loop contributions relative to the bare Coulomb term 
(i.e. the ratio of the second term to the first in the de- 
nominator of Eq. (3.22)). The RPA is generally expected 
to be valid when the loops dominate over other diagrams, 
i.e. ^ 1. Provided this condition is satisfied, we can 
also analyze the strong-coupling regime A ^ 1, and the 
crossover toward the weak-coupling one (A <C 1), i.e. we 
can hope to cover a wide range of a values. 

The calculated functions Fq and Fi in Eq. (3.23) are 



\/l~A2 



Fi(A) 



arccos A 



^+2A' 



2- A2 



In (a + v/A2~T) - 1 + 



i^o(A) 



A\/]~A2 
''-'.MX 



arccos A 



TT 



1 



A < 1, 

(3.25) 
A < 1, 

TT 

,A > 1. 



A\/A2~T \ ' / A 

(3.26) 

This leads to the system of RG equations for v and Z, to 
leading order in 1 /N 



dv 



dZ 



{F,iX)-Fo{X))v., 



FoiX)Z. 



At strong-coupling, A 1, one finds 



dv 



(3.27) 
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The first equation, after integration, leads to the low- 
energy result (fc — > 0) 



v(k)/v 



1] 



(3.31) 



which implies that the quasiparticlc dispersion is of the 
form 



w(fc) 



z = 1 - 



(3.32) 



The existence of the anomalous velocity dimension, 77, 
and consequently z 7^ 1, is characteristic of the strong- 
coupling regime Na — > 00 (Son, 2007). However this 
strongly-coupled fixed point is infrared unstable, since, 
due to the velocity increase, the RG for a fiows towards 
weak coupling. (One also expects that for certain N < 
Nc and a 3> 1 an excitonic gap can appear, which will 
be discussed in Section III.B.) In this regime Z can be 
approximated perturbatively (in 1 /N) as 



Z w 1 - -^^ ln{NaTT/i) ln(A/fc), 



Na > 1, (3.33) 



which can be obtained from Eq. (3.30) by ignoring the 
scale dependence of A. 

In the weak-coupling limit A 1, it is easy to verify 
that we recover the previous result (3.10) for the velocity 
V (leading to a flow for a towards zero), and the previ- 
ously encountered perturbative equation for Z 



dZ 
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' Ntt"^ 3 



1 



^a2ln(A/fc). (3.34) 



The last formula is written to first order in Na. 

Eqs. (3. 33), (3. 34) allow us to have a qualitative under- 
standing of the behavior of Z 0jS 0, function of the RG 
scale I. If the initial value of a is large, at the initial 
RG steps Z decreases logarithmically fairly fast (due to 
the weak ln(Q;) dependence in Eq. (3.33), even though a 
itself decreases). Eventually, when a has decreased sub- 
stantially {a ^ (ln(A/fc))~i), Z is governed by Eq. (3.34), 
meaning that Z will stop decreasing, and will level off for 
I = ln(A//s) 00. 

A numerical evaluation of the system of equations 
(3. 27), (3. 28) confirms the anticipated behavior and is 
shown in Fig. 11, (Gonzalez et at, 1999). (The equa- 
tion for the coupling A = ^Ne'^/{eov) is obtained by 
observing that {dX/dl) = {-l/v'^)^N{e'^ /eo){dv/dl), due 
to the lack of charge renormalization.) We conclude that 
the fiow of A is towards weak coupling, no matter how 
large its initial value is. Z does not renormalizc to zero at 
low energy due to the RG decrease of A. Thus, near the 
weak-coupling infrared fixed point, the marginal Fermi 
liquid (Eq. (3.19)) is ultimately not reached, and the 
system behaves as a Fermi liquid (although the quasi- 
particlc decay rate is non-Fermi liquid like, see below.) 
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FIG. 11 RG flow of the coupling A and the quasiparticle 
residue Z as a, function of the RG scale I; the infrared limit is 
at / — >■ 00. From Gonzalez et al, 1999. 



At higher energies however (away from the fixed point 
but still much lower than the bandwidth vA), the system 
exhibits marginal Fermi liquid behavior. 

At finite (but still small) density away from the 
Dirac point, i.e. fc 7^ 0, the logarithmic behavior in 
the infrared is cut-off by the Fermi momentum, i.e. 
ln(A/fci?), fcf/A — > 0, and the RG stops away from 
the fixed point. For comparison with experiments, the 
fiow toward this stable fixed point should be stopped at 
a scale set by the (small) density, temperature, or fre- 
quency, whichever is higher. 

One can also perform a numerical evaluation of the 
main RPA equation Eq. (3.21) (Pohni et al, 2007). For 
small density, and with logarithmic accuracy (ln(A/fc^)), 
this is equivalent to evaluating, by using the notation of 
Eq. (3.23), and taking into account Eqs. (3. 4), (3. 5), (3. 6) 



Z = {l-^T.^^'''^^/^u)- 



v* /v = Z [1 + 



8 



Fo{X)HA/kF) 

(3.35) 



Ntt^ 



FiiX)\n{A/kF) 



(3.36) 



Here v* is the renormalized velocity. At any finite den- 
sity the numerical evaluation of also picks up 
finite (sublcading) contributions, while it can be shown 
(Polini et ai, 2007) that the leading perturbative results 
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FIG. 12 (Color online) Exact evaluation of the RPA equations 
for (a) the quasiparticle residue, and (b) the Fermi velocity. 
On the horizontal axis / is defined as / = Na. A is in units 
of fci? . Values of A from ~ 10^ to 10^ correspond to density 
n from n ~ 10"cm"^ to n ~ lO^^cm"^ (while A ~ 10^ is 
ultra low density n ~ lO^cm"^). The values of A (in units 
of kp) can be converted into density n via: 



n = n/(10' 



The curves labeled 2DES refer to the 

and 



case of 2DEG with parabolic bands, where / = ^/2rs 
Ts ~ l/V^. From (Polini et ai, 2007). 



such as Eqs. (3. 33), (3. 34) are readily reproduced. The 
RPA results are shown in Fig. 12, and exhibit the natu- 
ral density dependence tendency, i.e. the strongest renor- 
malization occurs at the lowest densities. Similar RPA 
results have been obtained by Das Sarma et ai, 2007. 

A significant velocity enhancement was observed in the 
infrared conductivity (Li et ai, 2008), which reported 
around 15% increase of the Fermi velocity, having value 
as high as v* w 1.25 x 10®m/s at the lowest densities 
(compared to w sa 1.1 x lO^m/s at higher density). The 
system is at a finite Fermi energy fi sa 0.2eV. However 
the velocity renormalization is not logarithmic, and it is 
not clear what is the origin of this effect. 

A recent study of suspended graphene which measures 
the cyclotron mass (Ehas et al, 2011) has detected sig- 
nificant logarithmic renormalization of the Fermi veloc- 
ity, having the high value z;* « 3 x lO^m/s at the lowest 
densities n < lO^'^cm"^, almost three times the value 
at high density (n > 4 x lO^cm'^), Fig. 13(a). The 
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FIG. 13 (Color online) (a) Density dependence of the velocity 
for suspended graphene, from (Elias et al, 2011). The solid 
line is the result of KG treatment within RPA (Eq. (3.27)). 
(b) Reshaping of the Dirac cone due to the interaction-driven 
renormalization (increase) of the Fermi velocity at low mo- 
menta. The outer cone represents the linear Dirac spectrum 
without many-body effects. 



logarithmic renormalization of the velocity predicted by 
theory fits the data fairly well, and thus offers a direct 
proof that the Dirac cones can be reshaped by long-range 
electron-electron interactions near the Dirac point, as 
schematically shown in Fig. 13(b). Finally, ARPES mea- 
surements of quasi-frccstanding graphene grown on the 
carbon face of SiC have also detected logarithmic velocity 
renormalization (Sicgel et al, 2011). 



3. Quasiparticle lifetime 

The inverse quasiparticle lifetime (decay rate) due 
to electron-electron interactions, l/ree, is an important 
quantity which is relevant to many properties of graphene 
(and Fermi systems in general). In particular the depen- 
dence of l/Tee on energy (or temperature) determines 
the importance of the electron-electron interaction con- 
tribution, relative to other processes, to transport, and 
interpretation of spectroscopic features, such as ARPES. 

The decay rate is determined by the imaginary part of 
the self-energy, ImE(fc,w). The first diagram which has 
energy dependence, and thus a non-zero imaginary part, 
is the one bubble diagram of Fig. 7(b), whose real part 
is given by Eq. (3.16), i.e. behaves as in Eq. (3.19) at 
low energies. We can therefore deduce, for energies and 
momenta close to the mass shell (Gonzalez et ai, 1996), 

ImS(2'') (k, uj) - a^6{u} - vk) w, w « -yfc , (3.37) 

i.e. the decay rate is linear in energy. In addition, there 
is an on-shell ( "light cone" , w = vk) discontinuity, where 
the rate experiences a jump. This on-shell behavior is 
due to the fact that, for oj < vk, there is no phase space 
available for virtual interband particle-hole excitations 
(see Fig. 5), whereas such excitations arc possible for 
Lo > vk. 
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The above behavior is vahd at the Dirac point and T = 
fi ~ 0, while for small T, fi, it is valid for energies of order 
niax(T, /i). Notice also that the linear energy behavior of 
Eq. (3.37) is very different from the conventional Fermi 
liquid result ImS ^ (Das Sarma et al., 2007), which 
would occur for a finite Fermi surface (/i ^ 0) and is due 
to intra-band particle-hole excitations. 

The on-shcU discontinuity present at the one-loop 
level Eq. (3.37) disappears when the full RPA self- 
energy is evaluated (Fig. 10). In this case one obtains 
(Khveshchcnko, 2006) 

Iml]^'"^'^)(fc,w)-ln(7ra)6l(w-i)fc)(w-'yfc), cu f=i vk . 

(3.38) 

Away from the mass shell, the energy dependence is nat- 
urally linear: 

l^Y,'^RPA) ^) _ (^q,) cj > wfc . (3.39) 

The full dependence ImT,'- ^^^\k,uj) has to be evalu- 
ated numerically (Das Sarma et al., 2007), and the re- 
sults confirm the smooth rise of ImS'-^^'^) from the point 
Lo = vk. 

In the limit of zero doping /i — > 0, when the system ap- 
proaches the fixed point a = 0, we argued previously that 
the residue Z does not approach zero (i.e. the marginal 
Fermi liquid behavior ultimately docs not manifest itself.) 
On the other hand the marginal Fermi-liquid behavior is 
expected to be much more robust as far as the inverse 
lifetime, ImE ~ oj, is concerned, because the running of 
the coupling a(a;) only introduces logarithmic variation 
on top of a much stronger linear energy dependence. 

The linear decay rate discussed above is consis- 
tent with ARPES experiments (Bostwick et al., 2007; 
Zhou et ai, 2008), and STM measurements of graphcne 
on graphite (Li et ai, 2009a) (see also the discussion in 
Grushin et ai, 2009). 

B. Spontaneous mass generation 

It is an intriguing possibility that graphene can un- 
dergo a metal-insulator transition for strong enough 
Coulomb interaction a, due to an excitonic pairing mech- 
anism. We restrict ourselves to the charge neutrality 
point fi ~ since the excitonic pairing tendency de- 
creases quickly beyond that. 

1. Finite explicit mass 

Before we outline the main results, let us mention that 
an explicit gap can also open in graphene under certain 
conditions that depend on graphene's environment. For 
example there are suggestions of a detectable gap in situ- 
ations when graphene is on a substrate with specific sym- 
metry, creating sublattice asymmetry in the graphene 
plane, and thus making the graphene electrons massive 
(gapped) (Zhou et ai, 2007). Gaps can also be produced 



by confining the electrons into finite-size configurations, 
such as quantum dots (Ponomarenko et ai, 2008). In 
these cases the gap generation mechanism is not intrinsic 
to graphene, and the value of the gap depends strongly on 
the external conditions. However even in such situations 
interactions can play an important role by increasing the 
gap. 

Consider a gap arising from an external potential that 
alternates between the two sublattices 

n,nass = Ao ^ n^{R,) - Ao ^ n^{R^) . (3.40) 

Consequently an additional pseudospinor structure re- 
lated to 1T3 is generated, and the new Green's function 
has the form 

G{k,uj)^ rr^ vnr-v (^-^^^ 

uj — vcr ■ k — Ao(T3 — ij(fe, 1^) 

Here Aq is the explicit "mass" of the graphene elec- 
trons (while I](fe,w) contains the information about 
interactions, assumed to pcrturbatively renormalize 
all the other terms.) The new spectrum is then 
E{k) = ±^yv^k^ + Al, with a gap of 2Ao. Com- 
puting the Hartree-Fock interaction correction to Ag 
leads to a renormalized mass Aq (Kane and Mele, 2005; 
Kotov et ai, 2008a) 

Ao/Ao«l + |ln(i?/Ao). (3.42) 

The above enhancement can be substantial. For ex- 
ample for a bare gap due to spin-orbit coupling Aq ^ 
lO-^meV (Min et ai, 2006; Yao et al., 2007), and taking 
into account the bandwidth D = vA « 7eV, the logarith- 
mic factor is around 15. In fact one should integrate the 
RG equation for the renormalized mass Ag as a function 
of ln(A) simultaneously with the equation for the running 
coupling Q!(ln(A)), Eq. (3.13), down to the lowest infrared 
scale ^ Aq (bare gap). This leads to the stronger de- 
pendence Ao/Ao = (1 + f In (D/Ao))^ (Kane and Mele, 
2005), and the perturbative expansion of this result is 
Eq. (3.42). It is interesting to note that the logarithmic 
mass renormalization formula in graphene Eq. (3.42) is 
similar to the well-known expression for the electromag- 
netic mass of the electron (accounting for radiative cor- 
rections) in 3D relativistic QED (Weisskopf, 1939). 

2. Excitonic mass generation 

We now turn to the possibility of spontaneous gap gen- 
eration due to long-range Coulomb interactions (we set 
the explicit gap Ag = in Eq. (3.41)). In relativistic 
QED in two space (plus one time) dimensions, QED2+1, 
the study of this phenomenon, called chiral symmetry 
breaking, started quite a while ago (Appelquist et al., 
1986; Pisarski, 1984), and is still going strong today. 
Graphene is actually different from QED2+1 because only 
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FIG. 14 Schematic pliase diagram in the a — N plane. 

the fermions are confined to a 2D plane, while the field 
lines extend through the whole 3D space. In addition, the 
Coulomb interaction in graphene can be considered in- 
stantaneous since the speed of light c is much larger than 
the Fermi velocity (v w c/300). Hence, Lorenz invari- 
ance is not respected, which reflects the non-relativistic, 
purely band origin of the Dirac quasiparticles. The anal- 
ysis in relativistic QED reveals that dynamical mass can 
be generated below a critical number of fermion flavors 
Nc, with the mass scale set by the coupling itself, which 
has dimension of energy in pure QED2+1. A transition is 
also found in non-relativistic graphene, where the gener- 
ated mass scale is related to the ultraviolet energy cutoff 
(bandwidth D = vA) since the coupling a is dimension- 
less in this case. 

The gap equation can be obtained as a self-consistent 
solution for the self-energy within RPA (i.e. vertex 
corrections are neglected), and is referred to as the 
Schwinger-Dyson equation. It has the form 

/■ d^kduj V'^P'^ip ~k,e- uj)A{k, lo) 
^P'^)-'J (2^)3 uj^-v^k^-A^k)+tO+ ■ 

(3.43) 

The structure of the solution has been analyzed ex- 
tensively (Gamayun et al, 2010; Gorbar et ai, 2002; 
Khveshchcnko, 2009; Khveshchenko and Leal, 2004; 
Liu et ai, 2009) at different levels of approximation. 
The equation is simplified significantly if the static RPA 
potential is used l^^^'^(p, 0) (Khveshchenko and Leal, 
2004), while the dynamical equation has also been stud- 
ied on-shell (A(p, e = vp)) (Khveshchenko, 2009), as well 
as numerically (Liu et al., 2009). 

The mass gap A(p) has strong momentum dependence, 
due to the long-range nature of the Coulomb interaction. 
A(p) decreases at large momenta and reaches maximum 
value at small momenta where it levels off. For fixed 
physical value of iV = 4, a transition to a gapped state is 
found above a critical coupling ac- Some of the calculated 
values are: ac = 0.92 (Gamayun et ai, 2010), ac = 1.13 



(Khveshchenko, 2009). At strong coupling a —?' 00 the 
gap is non-zero only below a critical number of fermion 
flavors (since the effective interaction scales as 1/A^ in 
this limit); for example Nc « 7.2 (Khveshchenko, 2009), 
TVc « 7 (Liu et ai, 2009). 

Near the critical coupling the low-momentum gap 
scales as 

A(0) cx ZJexp f ^ - ) , (3.44) 

where C is a constant, the critical ae//,c = 1/2, and 
the form of the effective coupling aeff depends on the 
level of approximation used — for example an improve- 
ment over the static RPA potential leads to: aeff = 
a/{l + Nna/SVi) (which gives Nc ~ 7.2, a 3> 1, and 
ac = 1.13, iV = 4 (Khveshchenko, 2009)). The form of 
Eq. (3.44) suggests that the transition is of infinite or- 
der (Berezinskii-Kosterlitz-Thouless type). Even though 
Eq. (3.44) is only valid near the critical coupling, numer- 
ical results find that the gap in units of the bandwidth, 
A(0)/D, is exponentially small in a wide range of cou- 
plings (Khveshchenko and Leal, 2004). Since D « 7 eV, 
this implies A(0) ~ meV, i.e. a rather small gap value. 
Finally, recent work that takes into account the renor- 
malization of the coupling constant and the quasiparticle 
residue suggests that ac could be much larger than pre- 
viously found (Gonzalez, 2010; Sabio et ai, 2010a). 

The above results are based on various approxima- 
tion schemes and it is therefore important to com- 
pare them with direct numerical simulations of the lat- 
tice field theory model. Recent Monte Carlo calcula- 
tions (Drut and Lahde, 2009a,b,c; Hands and Strouthos, 
2008) provide strong evidence that spontaneous mass 
generation does occur, and give comparable values 
for the critical couplings: Nc ~ 9.6, a 3> 1 
(Hands and Strouthos, 2008), ac = 1.1, N = 4 
(Drut and Lahde, 2009b). Unfortunately the Monte 
Carlo simulations do not allow for an exact determina- 
tion of the gap size, and for that we can only rely on the 
previously described Schwinger-Dyson equation (leading 
to small gaps). For graphene deposited on SiO the value 
of asi02 ~ 0.79 and is therefore not enough to generate 
a gap; only experiments on ultrahigh mobility suspended 
samples can potentially reveal the insulating state. 

The overall phase diagram of graphene in the a — N 
plane is expected to look as shown in Fig. 14, with 
Q!c ~ 1 and Nc ~ 7 — 9. At finite temperature one 
expects the existence of a critical temperature Tc ^ 
A(0), while finite doping fj. very quickly destroys the 
gap (Liu et ai, 2009). Application of magnetic field 
perpendicular to the graphene layer leads to enhance- 
ment of the excitonic instability due to the formation of 
Landau levels (Gorbar et ai, 2002; Gusynin et ai, 2006; 
Khveshchenko, 2001a). In addition, it has been suggested 
that an in-plane magnetic field favors a gapped excitonic 
state (Aleiner et al., 2007), due to the instability of a 
system of electrons and holes polarized in opposite direc- 
tions. 
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FIG. 15 (Color online) ARPES data from (Bostwick et al, 
2010), showing strong features at the Dirac point, which is 
below the Fermi energy (at 0). The splitting shown in (H) is 
attributed to the presence of "plasmarons" — quasiparticles 
strongly bound to plasmons — and depends on the value of 
a (a ~ 0.5 fits the data.) 



The physical structure of the gapped state depends on 
the nature of pairing between the valleys — for example 
one can have charge density wave states (Khveshchenko, 
2001b) with modulation of the electronic density around 
the two sublatticcs (which corresponds to intravalley par- 
ing), or Kekule dimerization (Hou et al, 2007) which 
corresponds to tripling of the unit cell (intervalley pair- 
ing). One generally expects that interactions beyond 
the long-range Coulomb potential, such as short-range 
repulsion, would favor particular states, including time- 
reversal symmetry broken (spin) states. Further discus- 
sion appears in Section V.A. 



C. Finite density Fermi-liquid regime 

As the density increases above half-filling, i.e. 
graphenc is at a finite, not necessarily small, chemical po- 
tential /X, with a finite Fermi surface, a crossover towards 
a Fermi liquid regime takes place. In this case the lower 
(hole) band becomes irrelevant and the physics near the 
Fermi surface is dominated by intra-band transitions in 
the conduction (upper) band (assuming /i > 0) . However 
the physics near the Dirac point can still be very strongly 
affected due to the presence of plasmon and "plasmaron" 
features in the quasiparticle spectral function. 

The quasiparticle width near kp is quite simi- 
lar to the case of an ordinary 2D electron gas 
(Das Sarma et al., 2007; Hwang and Das Sarma, 2008b; 
Polini et al., 2008a), and is proportional to the second 



power of energy (or temperature), as in a Fermi liquid, 
while the quasiparticle residue is finite at the Fermi sur- 
face. 

The existence of a plasmon-related peak in the quasi- 
particle decay rate, which originates from intraband 
transitions in which an electron can decay into a plas- 
mon, was pointed out in the context of intercalated 
graphite, where the physics is dominated by graphene 
layers (Lin and Sliung, 1996; Shung, 1986b). For n-doped 
graphene (/i > 0), which is relevant to ARPES ex- 
periments, a double-feature is found in the decay rate 
ImS: a peak at positive energies, signaling an on- 
set of plasmon emission, and a sharp spectral feature 
at negative energies, below the Dirac point, and sepa- 
rated from it by an amount proportional to the plasmon 
frequency (Hwang and Das Sarma, 2008b; Polini et al, 
2008a). This is the so-called "plasmaron" — a resonance 
which consists of a quasiparticle strongly coupled to plas- 
mons (Lundqvist, 1967). Plasmaron features have been 
previously detected for example in optical measurements 
of Bismuth (Tediosi et al, 2007). 

The above calculations were done within RPA theory. 
Line widths have also been analyzed via ab-initio many- 
body methods (Park et al, 2009; Trevisanutto et al, 
2008). Experiments generally show a well-pronounced 
linear quasiparticle spectrum (Bostwick et al, 2007; 
Sprinkle et al, 2009; Zhou et al, 2007, 2008), with ad- 
ditional features near the Dirac point which seem to de- 
pend on the way graphene is prepared, and its purity. 
For example, gap-like features have been observed near 
the Dirac point (Zhou et al, 2007), and attributed to 
external, substrate-related factors. Bending of the Dirac 
spectrum (kink-like feature) was attributed to plasmons 
(Bostwick et al, 2007). Most recently manifestations of 
the sharp plasmaron spectral intensities have been ob- 
served in quasi- freestanding graphenc (Bostwick et al, 
2010), where a reconstruction of the Dirac point crossing 
seems to take place, as shown in Fig. 15. A diamond- 
like shape appears due to crossing of charge and plas- 
maron bands. Comparison of the RPA calculation for 
the energy splitting with experiment leads to the value of 
a K, 0.5 (Fig. 15.) Bostwick et al, 2010 also suggest that 
the plasmaron features were obscured in earlier measure- 
ments on non free-standing graphene (Bostwick et al, 
2007), due to the several times stronger screening (and 
consequently smaller a.) Perhaps most importantly, all 
the current activity in ARPES on different graphene sam- 
ples reveals that the electron-electron interactions can af- 
fect strongly the physics around the Dirac point, even for 
relatively large density (Fermi energy). 

Tunneling spectroscopy measurements, combined with 
ab-initio calculations, have also found evidence for 
density-dependent interactions effects in the tunneling 
current (Brar et al, 2010) which arise from the sharp 
spectral features in the quasiparticle decay rate below 
the Dirac point, as discussed above. 
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FIG. 16 (Color online) Inverse compressibility, measured by 
Martin et al., 2008. The red line is the compressibility of non- 
interacting Dirac fermions. 



D. Physical Observables 

The interaction-driven singular logarithmic structure 
near the Dirac point (for /Lt ~ 0) encountered in the 
fermion self-energy, and in particular the renormaliza- 
tion of the Fermi velocity, can manifest itself in numer- 
ous physical observables, such as the charge compressibil- 
ity and the spin susceptibility, which exhibit non Fermi- 
liquid behavior. Interactions can also affect the conduc- 
tivity near the Dirac point, leading to deviations from 
the celebrated quantized value co = e^/Ah expected for 
free Dirac fermions (Castro Neto et ai, 2009a). 



1. Charge and spin response 

a. Compressibility. First we discuss the compressibility 
K, which was recently measured (Martin et ai, 2008), 
Fig. 16, and it was concluded that no interaction effects 
were clearly visible in those samples. Theory predicts 
significant (a dependent) deviations from the free elec- 
tron behavior (Barlas et ai, 2007; Hwang et ai, 2007; 
Pohni et ai, 2008b; Sheehy and Schmalian, 2007). 

The computation of the compressibility requires knowl- 
edge of the ground state energy, which contains the first 
order Hartree-Fock exchange contribution E^x, and the 
correlation energy Ecorr, describing all the higher order 
effects. Keeping in mind applications of the theory for 
fairly strong coupling (a ~ 1), the contribution of Ecorr 
can be substantial. The correlation energy can be readily 
calculated within the RPA approximation, i.e. we take 
Ecorr ~ EjiPA- The total ground state energy E, per 
unit area, is the sum E = Ekm+Eex+EupA- The kinetic 
energy Ekm = {2/'i)vkpn, and n = (kpY /i: is the parti- 
cle density. The inverse compressibility is then calculated 
as 1/k = d^E/dn^, which is equivalent to the usual def- 
inition involving the variation of the chemical potential 
with density, (1/k) = d^/dn. For free Dirac particles 



this gives (I/kq) = W\/7r/(4n) — behavior which can be 
clearly seen in experiment Fig. 16. 

The interaction effects in the ground state energy ac- 
quire divergent contributions in the limit of small den- 
sity kp/Afu 0, similarly to the previously discussed self- 
energy (velocity) renormalization. Ignoring any finite 
(non-diverging) terms, one finds (Barlas et ai, 2007) 



Eex/n = -{vkp)hi{A/kp), (kp/A) ^ 0, 
6 



Erpa/ii = —G{a){vkp)ln{A/kp), 



6 



(3.45) 



(3.46) 



where the function G{a) is defined as G(a) ~ 
(1/2) dx{l + x'^)-^{Vx^ + l + NTra/8)-\ and, in par- 
ticular, at zero coupling G(0) = 1/3. The above results 
exactly follow the velocity renormalization, i.e. are equiv- 
alent to the substitution v — > v{kp) in the free compress- 
ibility (I/kq) = vy/ir/ (4n), where v{kp) is the running 
velocity calculated within RPA at the infrared scale kp. 
The result is particularly simple at the Hartree-Fock (ex- 
change) level (when the velocity follows Eq. (3.9)): 



- = v\ (I + -rHA/kp) + Oia' 
K V 4n V 4 



(3.47) 



and was obtained by a number of authors (Barlas et al, 
2007; Hwang et ai, 2007; Sheehy and Schmalian, 2007). 

The above results are valid at zero temperature. We 
also point out that exactly at zero density kp = 0, but 
T 7^ 0, the compressibility behaves as: ~ {v'^ /T){1 + 
(a/4) ln(ro/r))^, where Tq is the temperature related 
to the ultraviolet cutoff; since Av « 7 eV, then Tq « 
8 X lO^K. This is easily understood since in the infrared 
limit near the "critical point" n = T = it's the larger 
scale, either vkp, or the temperature T, which enters the 
physical observables (Sheehy and Schmalian, 2007). 

Of course Eqs. (3. 45), (3. 46) are valid only asymptoti- 
cally (kp 0), and at any finite density the compressibil- 
ity should be calculated numerically. This was achieved 
by expressing the ground state energy via the charge re- 
sponse function (Barlas et ai, 2007). 

Fig. 17, upper panel, illustrates the variation of 1/k 
with density for fixed interaction. Most notably, 1/k 
is larger than the free value 1/kq. Also, the full 
RPA implementation weakens the first order Hartree- 
Fock (exchange) result, due to the different signs in 
Eqs. (3.45), (3.46). For example, at a = 0.8 the RPA 
term is approximately 1/2 of the exchange, and thus has 
to be taken into account (although the RPA effects be- 
come weaker for a — > 0). Asymptotically, (k~^/kq ^) ^ 
ln(A/fci?), as kp/A —J- 0. The lower panel gives the vari- 
ation k/kq as a function of the interaction for different 
densities; naturally the deviation from the free limit in- 
creases with increasing interaction and decreasing den- 
sity. 

The increase of the inverse compressibility, kq/k, as 
a function of the interaction a (at fixed density), and 
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of the full RPA analysis led us to conclude that a < 0.1. 
It has also been argued that exchange and correlation ef- 
fects vanish and do not manifest themselves at all in the 
compressibility (Abergel et ai, 2009). These discrepan- 
cies indicate that the issue is still unsettled, while it's 
also possible (indeed, quite probable) that interaction ef- 
fects are obscured by charge inhoniogencitics (electron- 
hole puddles) in these samples. Nevertheless theory pre- 
dicts strong systemic (albeit logarithmic) deviations from 
Fermi-liquid theory, and it would be important to test 
these predictions in cleaner, more uniform, high-mobility, 
low-density samples. 
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FIG. 17 Upper panel: Inverse compressibility calculated at 
different levels of approximation as a function of density. The 
inset enlarges the low-density region. Lower panel (adapted 
from Barlas et al, 2007): Compressibility calculated within 
RPA, relative to tlie free level for different couplings and den- 
sities. Here A'' = 4 is the Dirac fermion degeneracy. The 
numbers refer to the values of A/kp, which can be converted 
into density n via: A/kr « 220/\/h, n = n/(10^°cm"^). This 
implies (A/kp) ~ 10^ for n ~ 10"cm"^ and (A/kp) ~ 10 for 
n ~ lO^cm"^ 



b. Spin susceptibility. The paramagnetic spin suscepti- 
bility, Xs , shows behavior very similar to the charge com- 
pressibility, i.e. {Xs/Xs,o) decreases as the interaction in- 
creases (Barlas et al., 2007). This is again related to the 
fact that x7^ is calculated via the ground state energy, 
and is proportional to the Fermi velocity v. It was also 
pointed out that the same effect, i.e. the logarithmic 
growth of the exchange energy, Eq. (3.45), can lead to 
suppression of ferromagnetism in graphene at low densi- 
ties (Peres et al., 2005). The full calculation of Xs within 
RPA was carried out by Barlas et al., 2007. 

On the other hand the orbital diamagnetic susceptibil- 
ity, Xdia, is proportional to v'^, because the quasiparti- 
cle current that couples to the vector potential contains 
V (the magnetic field is perpendicular to the graphene 
plane). Therefore interaction corrections lead to an in- 
crease of Xdia (Sheehy and Schmalian, 2007) and, conse- 
quently, orbital effects are expected to dominate in the 
susceptibility. At the Dirac point, hp = 0, one finds at 
finite temperature 

Xdia/Xd^a,0 = (l + f HTq /T)) ' , (3.48) 



with decreasing density (for fixed interaction) , represents 
non-Fermi liquid behavior, and reflects the lack of screen- 
ing. By contrast, in a 3D (and 2D) Fermi liquid with 
a screened potential kq/k decreases; for example within 
Hartree-Fock, kq/k w 1 — r^/G < 1, and eventually goes 
through zero, signaling an instability (Mahan, 2000) (al- 
though the critical value of 7's depends strongly on the 
level of approximation.) Such an instability does not oc- 
cur in graphene, which is related to the impossibility of 
Wigner crystallization (Dahal et al., 2006). It should be 
noted that for larger densities (larger than the density 
range shown in Fig. 17) the logarithmic corrections be- 
come unimportant and the system recovers the Fermi liq- 
uid behavior, i.e. eventually k/kq becomes larger than 1. 

Fits of the experimental data for k with adjusted 
(slightly larger) velocity v = l.l x lO^m/s show that 
a « (Fig. 16), while the use of w = lO^m/s by 
Sheehy and Schmalian. 2007 at the Hartrec-Fock level 
produced a sa 0.4. On the other hand, the application 



where the non-interacting Xdia.o = — e^w^/(67rc^T) 
(Ghosal et al., 2007). Here c is the speed of light. At T = 
0, 71 7^ 0, we have Xdia,o ^ —e^vl((?\/n), and interaction 
corrections readily follow from the v dependence. This re- 
sult is, strictly speaking, valid for T B ^ \i = v^/wn, 
whereas for _B = the orbital susceptibility is zero for 
/i 7^ as r — > 0, and is finite only when the Fermi energy 
is at the Dirac point. It has been suggested that an in- 
teraction driven positive (paramagnetic) contribution to 
the orbital susceptibility can therefore become dominant 
in doped graphene, Xorb ^ [e^v^/(/ic^)]a| lna|, a ^ 1 
(Principi et al, 2010). 



c. Specific heat. The specific heat is logarithmically sup- 
pressed due to the suppression of the DOS ~ v^^. Con- 
sequently Cv ~ Cv,o/(ln(To/T))2, T/Tq < 1, where 
Cv.o ^ T^/w^ is the free Dirac fermion specific heat. 
The full RPA calculation, valid also for large coupling, 
was carried out by Vafek. 2007. 
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d. Graphene as a quantum critical system. A uni- 
fied view of the above behavior is presented in 
Sheehy and Schmalian (2007), where it was stressed that 
the logarithmic corrections are manifestations of scal- 
ing behavior around the quantum critical point at n = 
0,T = 0. As discussed previously, at finite chemical 
potential, T = 0, n 7^ 0, graphene behaves as a Fermi 
liquid, whereas at T ^ 0, a quantum critical region 
fans out of the point n = 0,T = 0. In the critical re- 
gion it is natural to call graphene a Dirac liquid, where 
the proximity to the Dirac point is important for phys- 
ical phenomena at finite T. This puts graphene's be- 
havior into the general framework of quantum critical 
phenomena (Sachdev, 1999). In practical terms, it im- 
plies that the logarithmically divergent velocity contri- 
butions arc cut-off by the largest scale: temperature T, 
kp ^/n., or magnetic field. Computing physical quanti- 
ties in perturbation theory (Hartree-Fock or RPA) nat- 
urally involves these infrared scales. The separation be- 
tween the Dirac liquid and the Fermi liquid regimes in 
the n — T plane is defined by the crossover tempera- 
ture T*{n) = vkpil + (a/4) ln(A/fcF)), kp = yfwa, and 
thus the temperature dependencies quoted previously, 
are vahd for Tq > T > T*(n). The ultraviolet tem- 
perature scale To « 8 X lO^K, while for typical graphene 
densities n < IQi^cm-^ T*{n) - lO^K. 



2. Conductivity 

The behavior of the electrical conductivity in graphene 
has been extensively reviewed (Das Sarma et ai, 2011; 
Peres, 2010). It is believed that charged impurities 
and resonant scatterers are the main sources of scat- 
tering away from the Dirac point, and to extent the 
long- or short-range part of the Coulomb potential con- 
tributes to scattering is a matter of ongoing debate 
(Chen et ai, 2008; Monteverde et ai, 2010; Ni et at, 
2010; Ponomarenko et ai, 2009; Reed et ai, 2010). 

Here we will only mention effects related to long-range 
electron-electron interactions near the Dirac point. Inter- 
action corrections to the minimum metallic conductivity 
of free Dirac fermions, ctq = e^/{4h) = /h (Fradkin, 
1986; Lec, 1993), are more involved, because this expres- 
sion does not contain the quasiparticle velocity, while the 
electric charge is not renormalized. The debate was fu- 
eled in part by electrical measurements of the minimum 
conductivity (at the Dirac point) which turned out to 
be somewhat larger than ctq (Geim and Novosclov, 2007; 
Tan et ai, 2007). Theoretically, at T = (or T < a; 
where lo is the external frequency), it is expected that 
any interaction effect should have sub-leading character, 
and the frequency can enter only through the running 
of the coupling a{uj). Even though some debate still ex- 
ists (Herbut et at, 2008; Juricic et ai, 2010; Mishchenko, 
2008; Sheehy and Schmalian, 2009) as to the implemen- 
tation of the cut-off rcgularization procedure, the con- 



ductivity should have the form 

^('^)/^° = i + TTfRA^' ^'-''^ 

where the constant C w 0.01, as argued by Mishchenko, 
2008; Sheehy and Schmalian, 2009. The smallness of 
C reflects the near cancellation of self-energy and ver- 
tex corrections, and thus the effect of interactions is 
small. This value is also consistent with optical measure- 
ments on suspended samples (Nair et al., 2008), as well 
as graphene on a substrate (Li et ai, 2008), which find 
tT(a;) to be very close to CTq, and frequency independent 
in a wide range of energies. 

In the strict DC limit uj — 0, the presence of disorder, 
in combination with interactions, can alter the conduc- 
tivity. For example, for weak gauge field disorder (A) 
where an attractive line of fixed points exists (Fig. 9) 
with a* = ;|A, calculations show that the conductivity 
(on the fixed line) increases relatively to the free limit 
(Herbut et al., 2008): a = [n/2 + (4 - TT)A]e'^/h. For 
stronger scalar and vector disorder/interactions where 
the couplings run away to infinity the problem is non- 
perturbative, and a complex variety of behavior is ex- 
pected (Foster and Aleiner, 2008). 

For clean graphene at ^ = vkp = it was pointed out 
(Fritz et ai, 2008; Kashuba, 2008; MfiUer et ai, 2008) 
that at high temperature (compared to the frequency), 
the conductivity is expected to have the form: 

n 7fi 

a = ^^, Ta2»w, (3.50) 

where q:{T) = 4/ln(Aw/T) is the running Coulomb cou- 
pling. This form reflects electron-electron inelastic colli- 
sions with scattering rate l/ree cP'T . The linear tem- 
perature dependence is characteristic for Dirac particles. 
The above formula is valid as long as l/ree is the dom- 
inant scattering mechanism (collision-dominated trans- 
port), and implies that clean graphene at the neutral- 
ity point should exhibit a universal, interaction-limited 
conductivity, reflecting essentially the quantum critical 
behavior of graphene in this regime (T ^ [l). With 
increased doping {^/T), a crossover takes place to a 
Fermi liquid regime with screened interactions, where 
T-} - o?T'^l[i, (Miiller et al, 2008) and the conductiv- 
ity is dominated by charged impurity scattering. 

It has also been pointed out that for /i = graphene 
behaves as an almost "perfect" fluid, in a sense that 
its shear viscosity, 77, relative to the entropy density 
s is anomalously small: 77/5 — {0.13/a'^{T)){h/kB) 
(Miiller et ai, 2009). This ratio measures how strongly 
the excitations in a fluid interact. At room temperature 
77/ s of graphene is smaller than rj/s of any known cor- 
related quantum fluid, and is close to the lower bound 

IttF" proposed to exist for a large class of strongly 
interacting quantum field theories (Kovtun et al, 2005). 
Therefore, due to its quantum critical nature near the 
Dirac point, graphene is suggested to behave as a strongly 
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correlated quantum liquid and should exhibit signatures 
of electronic turbulence (Miiller et ai, 2009). 

E. Overview of main results 

Before we proceed with further topics related to inter- 
actions in graphene, let us broadly summarize the main 
findings and questions raised so far: 

1. For clean graphene at the neutrality point ^ = 0, 
interactions are not screened and are marginally 
irrelevant; the fixed point a* = is approached 
logarithmically (or, equivalently, the quasiparticle 
velocity increases logarithmically). From a theory 
standpoint, the approach towards this fixed point 
is well understood both from weak and strong- 
coupling (RPA) perspectives. Since in graphene 
one can have a ^ 1 under rather conventional ex- 
perimental conditions, our understanding of RPA 
calculations is important. RPA is justified only 
in the limit of large number of fermion species 
{N ^ 1), while for = 4 it should work for 
weak to moderate coupling; however there are indi- 
cations, coming mostly from two-loop calculations, 
that vertex corrections are numerically small, and 
thus RPA should work well. Disorder generally 
drives the system away from the clean fixed point, 
towards finite or even strong coupling, depending 
on disorder type. 

2. The resulting behavior near the Dirac point is that 
of a non-Fermi-liquid with a quasiparticle decay 
rate which is linear in energy, and decreasing quasi- 
particle residue. All physical characteristics related 
to the quasiparticle velocity (which increases loga- 
rithmically) are affected, and predicted to exhibit 
systemic, interaction dependent, deviations from 
their non-interacting values as the Dirac point is 
approached, either as a function of density or tem- 
perature. 

3. Can graphene be driven into an excitonic insulating 
state? At the Dirac point the long-range Coulomb 
interactions can lead to bound electron- hole pairs, 
creating a gap. There has been intense debate 
whether this can happen under realistic conditions 
— since the critical interaction strength appears to 
be ac ~ 1, it seems possible to occur in suspended 
samples (a = 2.2). So far no experimental indica- 
tions have been observed. 

4. What is the value of the interaction a? Clearly, 
since a = 2.2/eo is dielectric constant depen- 
dent, working with different substrates could 
lead to changes in interaction-dependent effects 
(Jang et al., 2008). There are also suggestions that 
graphene has an "intrinsic" value of a (Reed et al, 
2010), arising from dynamical dielectric screen- 
ing. The polarizability of the Dirac fermions was 



found to be amplified by excitonic effects, improv- 
ing screening of interactions between quasiparti- 
cles. This analysis leads to values of a ranging 
from a sa 1/7 in the static limit to a ~ 2 at high 
frequencies. Very recent measurements of the cy- 
clotron mass in suspended graphene (Elias et al., 
2011) have found logarithmic velocity renormaliza- 
tion and extract, within the RPA scheme, an ef- 
fective value of graphene's dielectric constant ec ~ 
3.5. One can also expect that near the Dirac point, 
where interactions lead to singular effects, addi- 
tional factors can be important such as disorder, 
inhomogeneities, rippling, etc., and thus obscure 
the clean behavior. 

5. In the Fermi-liquid regime, where interactions are 
screened, the physics near the Dirac point can still 
be strongly affected — this is due to resonant fea- 
tures in the quasiparticle self-energy, reflecting in- 
teractions of quasiparticles with plasmons. 

IV. THE COULOMB PROBLEM AND CHARGED 
IMPURITIES 

The consideration of non-interacting Dirac electrons 
in 2D under a Coulomb field is of paramount relevance 
for graphene, and for several reasons. First of all, the 
Coulomb problem for relativistic fermions has many fea- 
tures that are unfamiliar in condensed matter systems, 
and which resemble long standing predictions made in the 
context of QED in strong fields. As such, and given that 
having a ^ 1 makes graphene intrinsically strongly cou- 
pled, it can provide the first experimental ground for test- 
ing many elusive predictions from strong-coupling QED. 

On the other hand, the single particle Coulomb prob- 
lem constitutes the first step in addressing nontrivial fea- 
tures of the full, many-body interacting problem. Char- 
acteristics like non-linear screening, or the supercritical 
instabilities, provide valuable insight in grasping some 
proposed many-body effects, like exciton condensation, 
or spontaneous mass generation in graphene. 

Historically, however, the motivation for studying the 
Coulomb problem comes from the seminal experimental 
observations (Novoselov et al. , 2004a) that the field effect 
in graphene prepared on Si02 is characterized by carrier 
mobilities that do not depend on the Fermi energy or 
carrier density (the DC conductivity, a = me|n|, with 
TO ~ const.), and that carriers are chiral Dirac fermions 
in 2D (Novoselov et al., 2005; Zhang et al, 2005). Early 
semiclassical investigations (Adam et al., 2007; Ando, 
2006; Nomura and MacDonald, 2007, 2006) showed that 
such linear-in-dcnsity conductivity could be explained by 
scattering of unscreened Coulomb impurities, which are 
typically seen in silica in concentrations of ^ lO^'^cm"^ 
(Ando et al., 1982). As a result, transport in the pres- 
ence of charged impurities rapidly became one of the 
most studied topics in the quest for the ultimate mobil- 
ity in graphene. Since, as we saw before. Coulomb's law 
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is exactly preserved in undopcd graphene, and approx- 
imately preserved for small and moderate doping, the 
scattering processes are essentially governed by the bare 
Coulomb problem, unlike conventional metals, where 
screening is perfect. A thorough understanding of this 
problem is therefore important not only for its theoreti- 
cal relevance and its import on electron-electron interac- 
tions, but also for its experimental implications, and our 
understanding of transport in graphene. 

Finally, it is highly significant that this is an exactly 
solvable problem. This means that most quantities can 
be obtained exactly, allowing us to unveil many interact- 
ing and non-interacting effects that are not within reach 
of the pcrturbative approaches already discussed. We 
proceed to show several such features. On account of 
the long range nature of the Coulomb field, inter-valley 
processes are not relevant, and hence we will solve the 
problem within each (independent) valley in the Dirac 
description of fermions in graphene. 



A. Exact Solution of the Coulomb Problem 

1. Wave Equations and Spectrum 

A Coulomb center of charge Z\e\ generates the poten- 
tial U{r) = Ze'^/{eor) for the electrons. Without any 
loss of generality let us consider Z > 0. The electronic 
dynamics is governed by the wave equation 



r 



cjsMv *(r) = £'^'(r) 



(4.1) 



Here we use g = Za = Ze^jie^v). with cq reflect- 
ing the effective dielectric constant of the embedding 
medium, and the mass M accounts for the more gen- 
eral possibility of a symmetry breaking gap. Through- 
out this chapter we shall use the scaled energy and 



mass £ = -B/i', m = Mw, and k = \J ~ rn^ . 
Even though m = for ideal graphene without in- 
teractions, nonzero m can be induced in many ways. 
One of them is through interaction with suitable sub- 
strates, of which some experimental hints have been 
reported (Griineis and Vyalikh, 2008; Li et at, 2009a; 
Martinazzo et ai, 2010; Zhou et ai, 2007). In terms of 
the original tight-binding Hamiltonian, the mass M aris- 
ing from a sublattice symmetry is related to the parame- 
ter Ao introduced in eq. (3.40) via Mv'^ = Aq. The axial 
symmetry of the potential allows us to use the eigen- 
states of the total pseudo angular momentum, = Lz + 
CTz/2, which is conserved (DiVincenzo and Mele, 1984). 



We write 



-l/2r 



, ^[F,(r)$,_i/2(0),zG,(r)$^-+i/2(</))], 
where j = ±1/2, ±3/2, . . . are the eigenvalues of J^, 
and the cylindrical harmonics read ^p{4>) = e^^'^/^/2Tr. 
A detailed derivation of the 2D Dirac equation for gen- 
eral radial potentials is given by Novikov, 2007a. In our 
case, Eq. (4.1) reduces to the following radial equations 



(Khalilov and Ho, 1998; Novikov, 2007a) 

[m-e- g/r] F, (r) + [9, + ]/r]G, (r) = (4.2a) 
[dr - j/r] Fj {r) + [m + e + g/r] Gj (r) = 0. (4.2b) 

This coupled pair of first order equations can be straight- 
forwardly reduced to two decoupled second order equa- 
tions. Free solutions [g — 0) of (4.1) exist when |e| > |m|, 
and are simple spherical waves whose fc-normalized ver- 
sion reads 




y'\ e + m\ J j-i/2{kr) $j-i/2 
ise\/\e - m\ Jj+i/2(kr) ^j+1/2. 



(4.3) 



(sj. = sgn(a::)). For nonzero g, one readily sees from (4.2) 
that the solutions at r ~ behave as 



F{r),G{r) - r 



±7 



7 



(4.4) 



The general exact solution is given in terms of confluent 
hypergeometric, or Whittaker's functions, both in the 
massive (Gamayun et ai, 2009; Gupta and Sen, 2008; 
Gupta et ai, 2010; Khahlov and Ho, 1998; Novikov, 
2007a; Pereira et ai, 2008a), and massless cases 
(Gupta and Sen, 2009; Pereira et ai, 2007; Shytov et ai, 
2007b). In the massless case, one can map (4.2) into 
the familiar Coulomb radial Schrodinger equation in 3D 
(Pereira et ai, 2007): 

d^rf± + [e' + 2ge/r - 7(7 t l)/r^] f±{r) = 0, (4.5) 

where the f± arc linear combinations of F and G, 
takes the place of the Schrodinger energy, and 7 plays 
the role of angular momentum. Since the solution is for- 
mally the same, the appearance of instead of e means 
that the massless case admits no bound solutions, as we 
expect on account of the absence of a spectral (mass) 
gap. The massive case, however, has a well defined infi- 
nite spectrum of bound solutions when \e\ < |m|, given 
by (Khalilov and Ho, 1998) 



(4.6) 



lowest level is given by Sg = £0,1/2 = Sgmi/l — (2g)^. 



2. Supercritical Instabilities 

Consideration of eq. (4.4) immediately reveals a com- 
plication a g > gc = 1/2; because 7 becomes imaginary 
for the lowest angular momentum channels {j = ±1/2). 
The solution (4.4) is neither regular nor divergent, but 
rather oscillates endlessly towards r = 0. This is patho- 
logical because the space of solutions is of dimension 2, 
and we can no longer discard an irregular contribution 
since both linearly independent solutions are square in- 
tegrable. In other words, there is no boundary condition 
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FIG. 18 (Color online) Schematic drawing of the level diving 
process in the supercritical regime, and of the resulting quasi- 
spectrum of levels for massive and massless fermions. 



at the origin to univocally select the solution. Secondly, 
in the massive case the level eq becomes imaginary, sig- 
naling a loss of self-adjointness of the Dirac Hamiltonian 
for g > 1/2. 

Physically, both effects are a symptom that the po- 
tential has such a strong divergence that particles arc 
inexorably attracted and "fall" into the origin, leading 
to a collapse of the system (for example, the endless os- 
cillations can be read as an infinite phase shift). This 
"fall to the center" is a general characteristic of diverg- 
ing potentials in any dimension of space. For power law 
potentials, one particular power signals the threshold of 
criticality. The Coulomb potential is the marginal case 
for the Dirac equation (both in 2D and 3D), just like the 
potential 1/r^ is the marginal case of the 3D Schrodinger 
equation (Landau and Lifshitz, 1981). This, of course, 
begs the question of regularization. Regularizing the po- 
tential introduces an additional boundary condition at 
some short distance R, which allows a formal solution, 
and cures the total collapse of the system (Case, 1960; 
Perelomov and Popov, 1970). In graphene the lattice is 
the natural regulator and there arc no ultraviolet issues. 
But the physics in the supercritical regime depends ex- 
plicitly on the short range details. 

This supercritical collapse has a long history in the 
context of QED, where the Dirac equation stands as the 
basis for understanding the stability of matter. In QED 
the collapse would occur for Zuqeb > 1, which lead to 
extensive investigations regarding the stability of heavy 
nuclei having Z > Z^ — 137 (Case, 1960; Greiner et at, 
1985; Popov, 1971a,b; Zeldovich and Popov, 1972). Af- 
ter regularization Zc — >■ 170, which makes the problem 
highly academic, and QED's predictions untestable. In 
graphene, on the contrary, Z^ ^ I, which opens the real 
possibility of testing the supercritical instability in a con- 
densed matter setting. 



a. Massive Electrons. To understand the physics in the 
supercritical regime we can follow the level as 
the coupling increases (Fig. 18) (Greiner et ai, 1985; 
Pereira et ai, 2008a; Zeldovich and Popov, 1972). For 



the pure Coulomb case, £0(9) decreases towards zero 
in a singular way at g = Qc- In a regularized poten- 
tial, £g depends also on the cutoff radius R, and is al- 
lowed to monotonically penetrate the negative energy re- 
gion, until eventually touching the lower continuum at 
e = —777. If g is further increased, Eg dives into the hole 
(positron) continuum and becomes a resonance. Other 
levels will sequentially follow at higher g. The diving 
point for scig) defines a renormalized critical coupling, 
gc > gc that is characterized by a log singularity at 
mR ^0: (jc — gc + T^^/iog^{mR) (Gamayun et al, 2009; 
Khalilov and Ho, 1998; Pereira et ai, 2008a; Zhu et al., 
2009), strongly depending on the regularization. 

This diving of bound levels entails a complete re- 
structuring of the vacuum. If the level was empty, an 
electron-hole pair will be immediately created: the elec- 
tron remains tightly bound and shielding the center, 
while the hole is ejected to infinity (Greiner et ai, 1985; 
Zeldovich and Popov, 1972). The supercritical regime 
is thus characterized by spontaneous pair creation, or 
a spontaneous Schwinger mechanism (Schwinger, 1951). 
One expected consequence is a strong signature of these 
resonances in the hole sector of the scattering and trans- 
port cross sections. 

An essential detail is that these resonances are not 
usual bound levels diluted inside a continuum, where 
their lifetime essentially disappears. One consequence 
of the chiral nature of Dirac fermions, combined with 
the long range tail of the Coulomb potential, is that the 
supercritical levels in the relativistic Coulomb remain 
sharply defined, with diverging lifetime. For example, 
for S states (j = 1/2), one shows that these resonances 
follow (Gamayun et ai, 2009) 
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when g ^ (jc, and where 



77, Pc = ^J~gl - 1/4. In 



real space the localization of the supercritical levels is 
controlled by the reduced Compton wavelength: Ac = 
1/(7777;). The modulus squared of their wavefunction de- 
cays as ^I'^^E' cx exp(— y^gr/Ac) and, consequently, even 
inside the continuum, such levels retain a highly localized 
nature, which is why they are so relevant, in particular 
in their potential for screening (Pereira et ai, 2008a). 



b. Massless Electrons. The spectrum in this case is con- 
tinuous everywhere, and thus there is no sequential div- 
ing and restructuring of the hole continuum as described 
above. But the pathology associated with Eq. (4.4) still 
exists. Physically, the massless situation is rather more 
catastrophic since the solution in a regularized potential 
reveals an infinite number of quasi-localized resonances 
in the hole sector (Gamayun et ai, 2009; Pereira et ai, 
2007; Shytov et ai, 2007b). This is a highly non-trivial 
effect for several reasons: (i) in the massless case there 
is no natural length scale in the problem to characterize 
such localized states; (ii) the system abruptly develops 
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an infinite quasi-bound spectrum at g > gc, when its 
spectral fingerprint is rather featureless for g < gc', (hi) 
the infinite spectrum has the potential to over-screen the 
Coulomb center. In addition, unlike the massive case, 
here the critical coupling remains unchanged at = 1/2, 
and no qualitative features (like how many, if any, states 
have dived) depend on the magnitude of the rcgulariza- 
tion distance. The spectrum of supercritical resonances 
behaves as (Gamayun et ai, 2009; Gupta and Sen, 2009; 
Shytov et ai, 2007b) 



{a,h)^0{g), (4.8) 



which has an essential singularity at g^ an energy 
scale/lower bound set explicitly by the regularization dis- 
tance, i?, and diverging lifetimes close to the critical 
point. Since the width of these states vanishes linearly, 
they are practically bound states (hence the designation 
quasi-bound states). In real space, the localization scale 
is determined by the regularization distance R itself. 

Since meso and nanoscopic devices arc of high interest, 
it is pertinent pointing out that massless Dirac fcrmions 
in a finite-sized system mimic in all aspects the physics 
of massive electrons, as a result of the linearly vanishing 
DOS and the effective gap coming from finite-size quan- 
tization (Pereira et ai, 2008a). 



3. DOS, Scattering and Transport Cross Sections 

Here and in the coming sections we shall be concerned 
mostly with massless Dirac fermions, except when explic- 
itly stated otherwise. The local density of states (LDOS) 
and cross sections are useful quantities insofar as they are 
directly accessible in local probe and transport experi- 
ments. The LDOS per unit area and spin is isotropic, and 
can be written in closed form in terms of partial waves 
as N{e,r) = J2j ^ji^,''^), (Pereira et ai, 2007) with 



nj{e,r)-- 
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(4.9) 



for g < gc, and Fi represents the Coulomb func- 
tion Fi{—gse, \e\r) (Abramowitz and Stegun, 1964). The 
function N{e,r) is plotted in Fig. 19(a) for different cou- 
plings and distances. Apart from the evident particle- 
hole asymmetry, the LDOS remains rather featureless, 
even at the shortest distances. It g > gc the corre- 
sponding analytical expression obtained in the regular- 
ized potential is more complex, but still has a closed form 
(Pereira et ai, 2007). In this case, supercritical channels 
(|j| < 1/2) need to be isolated from undercritical ones 
(|j| > 1/2). yielding two contributions to the LDOS: 

7V(£,r)= ^ n,{e,r)+ ^ n,{e,r). (4.10) 
b\<\a\ \3\>\9\ 

The total LDOS for this case is shown in Fig. 19(c) for 
g = 1.0, and at different distances to the impurity. It is 
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FIG. 19 (Color online) (a) LDOS, iV(e, r) at r = a for several 
couplings g < gc- The inset shows N{e,r) for g = 0.27 and 
different r. For comparison, the exact LDOS calculated in the 
full tight-binding lattice for the same parameters is shown eis 
dashed lines. In the horizontal axis the energy is in units of the 
hopping t. (b) The weak coupling transport cross section as 
a function of g. The inset shows the phase-shifts for different 
j. (c) LDOS, N{e,r) at several distances r, for g = 1 > 
gc- The inset shows the oscillating LDOS correction for e > 
0. (d) Energy dependence of the phase shifts (top) and the 
supercritical contribution nj{e,r) to the LDOS (bottom) for 
3= 1.0. 



now clear that strong resonances, decaying rapidly with 
distance, appear in the vicinity of the Dirac point, signal- 
ing the presence of the quasi-bound levels (Pereira et ai, 
2007; Shytov et ai, 2007b). Their exponential accumu- 
lation at £ = is confirmed in Fig. 19(d) where we show 
the supercritical contribution nj{e,r) as a function of 
log(|e|). At positive energies the LDOS exhibits peri- 
odically decaying oscillations in er [inset of Fig. 19(c)], 
with cxtrema separated by w nn, within logarithmic ac- 
curacy (Shytov et ai, 2007a). When directly measured 
in STM such oscillations can be used to extract the elec- 
tronic dispersion, as done by Ouyang et ai, 2002. 

We point out that, since the solution of the supercrit- 
ical problem involves a nontrivial ad-hoc regularization, 
these results have been checked numerically against ex- 
act solution of the full tight-binding problem in the hon- 
eycomb lattice, being found that the analytical Dirac re- 
sults reproduce the full lattice problem down to distances 
as small as the lattice scale (Pereira et ai, 2007). 

The striking differences between the two regimes and 
the violent modification of the ground state at strong cou- 
pling arc likewise evident in the behavior of the scattering 
phase-shits, Sj{e). They admit closed formed expressions 
at both g < gc (Novikov, 2007a; Pereira et ai, 2007; 
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Shytov et al, 2007a) and g > gc (Castro Neto et ai, 
2009b; Shytov et al, 2007b). For example, the under- 
critical S'-matrix reads (Novikov, 2007a) 

S,ie)-^-"<-> ^^;""\ . (4.11) 

7 — tgs^ i (1 + 7 + tgse) 

which is energy independent, but considerably asymmet- 
ric with respect to the sign of g. The corresponding Sj 
are shown in Fig. 19(b) (inset) as a function of coupling 
strength. Note how S1/2 (the most important partial 
wave) behaves rather differently from the others: only 
61/2 shows the expected sign for the attractive/repulsive 
situations. On the other hand, in the supercritical regime 
there is a strong e-depcndence of Sj . In the top row of 
Fig. 19(d) we present {6j mod tt) as a function of log(e). 
In the attractive sector (e < if 5 > 0) the abrupt steps 
centered around tt/2 mark the position of the infinite 
quasi-bound spectrum (which, as per (4.8), accumulates 
exponentially at e = 0), whereas in the attractive sector 
6j{e) is smooth. 

Knowledge of the phase-shifts allows direct calcula- 
tion of the full transport cross-sections for our 2D Dirac 
fermions: 

Atr(e) - ^^sin2(5,+i/2(e) -5,_i/2(£)) (4.12) 



(Katsnelson, 2006; Novikov, 2007a). The profile of At,, x e 
at weak coupling is shown in Fig. 19(b). When scat- 
tering is due only to unscreened charges, the marked 
asymmetry between g > and g < can be used to 
extract the density of positively and negatively charged 
impurities (nf) from a single measurement of the elec- 
trical conductivity, cr, as a function of carrier density 
(Novikov, 2007b). This technique has been used in some 
experiments (Chen et ai, 2009a,b, 2008), but the asym- 
metry effect can be easily masked by other spurious in- 
fluences (Barraza-Lopez et ai, 2010; Huard et ai, 2008; 
Nouchi and Tanigaki, 2010). Moreover, on account of the 
e-independence of 6j in (4.11), the corresponding Drude 
conductivity, a = 47re^/i/(uniAtr/i^), is immediately seen 
to scale linearly with density: a (X fj,'^ cx n. Therefore, 
the linear-in-density conductivity, which appears already 
in the first Born approximation, remains when the cross 
section is calculated exactly. 

For supercritical potentials, and similarly to the 
LDOS, there will be undercritical and supercritical par- 
tial waves contributing to Atr(e) [cfr. eq. (4.10)]. The 
latter give rise to strong peaks in the transport cross- 
section at densities for which the Fermi energy matches 
the levels £„ (Shytov et al. , 2007b) , tallying with the be- 
havior of the DOS. 



B. Induced Charge and Screening 

First attempts at understanding screening in graphene 
date back to DiVincenzo and Mele, 1984, where it was 



recognized that conventional procedures of the theory of 
metals, like self-consistent screening, linear response or 
Friedel sum rules, are not straightforward in this sys- 
tem. For example, within the Dirac (effective mass) ap- 
proximation, the ultraviolet cutoff scale enters explic- 
itly in Friedel's sum rule, and Levinson's theorem is 
modified (Lin, 2006) (Levinson's theorem is one of the 
fundamental results in quantum scattering theory, as- 
serting that in the Schrodinger's equation with a non- 
singular spherically symmetric potential the zero en- 
ergy scattering phase-shift exactly counts the number of 
bound states: Si{0) = Nin). One consequence is that 
a naive application of Friedel's sum rule can yield di- 
vergent displaced charges (DiVincenzo and Mele, 1984). 
Even though these divergences are artificial in the target 
lattice problem, they point, already at a single particle 
level, to the anomalous screening properties of graphene. 



1. Weak Coupling [g < gc) 

a. Non-interacting Induced Charge. Knowledge of the ex- 
act LDOS within the Dirac approximation (Sec. IV. A. 3) 
allows the straightforward calculation of the perturbation 
to the electronic density induced by the Coulomb center. 
The induced density is defined as 6n{r) = n{r) — n°(r), 
and is related to the LDOS via (for undoped graphene at 

zero temperature) n{r) = ^j(^) = 'J2j I-d "jC^' ^)ds, 
where D is the cutoff scale for the linearly dispers- 
ing band. The induced charge density is just Sp{r) ~ 
-'\e\5n{r). Closed form expressions for nj{r) are pro- 
vided in (4.9). One difficulty with this approach is that 
the resulting density per partial wave behaves asymptot- 
ically as 



5nAr — )• 00)^ - 
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(4.13) 



which diverges upon summation over j (a reminiscence 
of the problems associated with the ultraviolet scale al- 
luded to above). In the above expression D and 
represent the cutoff in the presence and in the absence 
of the coulomb center, respectively. Since the sublead- 
ing terms in (4.13) are convergent in ), we regularize it 
by taking a position dependent cutoff: D — >■ + ^. 
As a result, the total induced density acquires the form 
5n{r) ^ H{D^r)/r^ , where H{x) is a constant-amplitude 
oscillating function (Pereira et ai, 2007). Since it is de- 
sirable to have control over the validity of the regular- 
ization procedure outlined above, we have calculated the 
total induced density 5n{r) in the full tight-binding prob- 
lem, via exact diagonalization. The result is plotted in 
Fig. 20(a), and unequivocally shows the predicted l/r^ 
decay, with oscillations on the scale of the lattice. Such 
fast decay implies that the induced charge concentrates 
within a small vicinity of the impurity. Moreover, the 
numerical results in the lattice further suggest that such 
distance is of the order of the lattice parameter a: the 
inset in Fig. 20(a) reveals that the total charge pulled in- 
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FIG. 20 (Color online) (a) Induced electron density, iS?i(r), 
plotted as a function of distance to the Coulomb center, for 
different impurity strengths, g < g^. Data obtained from full 
diagonalization of the tight-binding Hamiltonian in a lattice 
with 124^ atoms. Black lines are oc 1/r'^, and guides for the 
eye. Inset shows the saturation of the integrated charge ac- 
cumulated inside r < i?max, as a function of -Rmax- (b) Same 
as (a), but for the supercritical case, g > gc, and the dashed 



line is now oc 1/r 



side a region r < i?max saturates within very few lattice 
spacings. In fact, since D° oc 1/a, in the limit a — > 
(where the effective mass description is meaningful) the 
analytical expression Sn{r) ~ H{D'^r)/r^ can be seen as 
a representation of the 2D Dirac-delta function. In other 
words, we expect the induced charge density to behave 
as 



dp{r) ~ — |e| Sn{r) 



a-5-O 



-Q\e\Sir). 



(4.14) 



The same conclusion follows from a modified Friedel ar- 
gument (Shytov et ai, 2007a); and from the exact cal- 
culation of the non-interacting Green's function in the 
Coulomb field (see below) (Terekhov et ai, 2008). The 
induced charge has a screening sign, as expected, but 
the strongly localized distribution of the induced charge 
(4.14) implies that undoped graphene cannot screen 
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FIG. 21 (Color online) (a) Total integrated charge in the 
vicinity of the impurity, Q, obtained: from exact diagonal- 
ization in the lattice (dots), from RPA (4.16) (blue), and 
from the exact Green's function in the Coulomb field (4.19) 
(red), (b) The self-consistent Zed, obtained from Eq. (4.21) 
(Terekhov et al., 2008). Numerical data (dots) is plotted after 
accounting for finite-size renormalization of gc (Pereira et ai, 
2008a). 



in the usual sense, because it merely renormalizes the 
strength of the impurity: Z — ^ ZcS = Z — Q. This 
leaves Coulomb's law unaltered, except for the substitu- 
tion Z — > Zcff. 



b. Linear (RPA) Screening. Single particle results, like 
the one above, are not generally sufficient to draw con- 
clusions about screening. Consider now the same prob- 
lem in linear response, at the RPA level, which is jus- 
tified for small, undercritical couplings. Within the 
RPA, the Fourier transform of the statically screened 
potential is given by Us{q) = Uo{q)/[l - n^'^\q)V{q)] 
(Fetter and Walecka, 1971), where V[q) = 2T:e^/{eoq) is 
the electron-electron interaction, and UQ{q) ~ ZV{q) the 
external impurity potential. From (2.14) we know that 
n(^)(gf 0) w -g/(4w), and hence 



Us{q)-Uoiq){l + 



Uo{q) 
erpa 



(4.15) 



Therefore linear response confirms the absence of screen- 
ing, except for the trivial renormalization of the static 
dielectric constant: eg — >■ erpa = eo(l + Tra/2) (Ando, 
2006). Likewise, the induced density can be computed 
in linear response from dn{q) = — ZV {q)Il{q) or, in the 
RPA: 



Sn{r) 



Z / dq 



l-n(i)(g)F(g)' 



(4.16) 



yielding Sp{r S> a) ~ -~5{r)Z\e\TTa/2 to hnear order 
in a (Kolezhuk et ai, 2006). This is exactly what was 
obtained in (4.14) from a single particle, wavefunction, 
perspective. In addition, the argument that the Fourier 
transform of dn{r) is dimensionless can be used to show 
that it should be a pure constant in undoped graphene, 
for which there is no natural length scale. As a result. 
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that Sp{r) cx S{r) remains true in all orders of pertur- 
bation theory (Biswas et ai, 2007). For consistency, the 
total induced charge, Q, introduced in (4.14) is then given 
by 



(higher orders in Za). 



(4.17) 



To verify this correspondence we can compare (4.17) with 
the value of Q extracted from the non-interacting exact 
diagonalization in the honeycomb lattice. As shown in 
Fig. 21(a), the numerical Q for different values of Z fol- 
lows the relation (4.17) for most of the range < g < gc, 
thereby confirming the correspondence, and showing how 
weakly undoped graphene screens (Pereira et al., 2007; 
Shytov et ai, 2007a). Given that only the global di- 
electric constant is affected, one can say that undoped 
graphene screens like an insulator. 

At finite densities, however, the system screens like 
a conventional metal. This derives at once from the 
fact that, at finite Fermi momentum, Il^^^q « 0) « 
— 2fcF/(7rw), no longer vanishing, and leading to the 
screened potential 
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(4.18) 



qs = Aakp playing here the role of inverse screen- 
ing length (Ando, 2006; Nomura and MacDonald, 2006). 
Contributions from interband transitions can be simply 
incorporated by renormalizing the background dielectric 
constant by the factor (1 -I- 7ra/2), as in eq. (4.15). Us- 
ing (4.16), the total integrated charge is now seen to be 
J 5p{r)dr = —Z\e\. This means that, unlike the un- 
doped situation, at finite electron densities the system 
completely screens the Coulomb center, just as expected 
in a metallic system (Castro Neto et al., 2009b). 

For transport considerations it is important to under- 
line that, even though at finite densities charged impu- 
rities have a finite range determined by qs, the Boltz- 
mann conductivity remains linear in density. This hap- 
pens because the screened potential (4.18) entering in 
the relaxation time calculation, maintains the same de- 
pendence with kp. From this perspective, the mobility 
remains constant in density for both screened and un- 
screened charges, differing only by an overall constant 
related to eRPA(fcF) (Nomura and MacDonald, 2006). 



c. Nonlinear Screening. As Fig. 21(a) documents, even as 
linear response is acceptable at small values oi g = Za, 
the approximation becomes increasingly unwarranted as 
g nears the critical threshold, gc = 1/2, which is non- 
perturbative. Rather than analyze this limit on the ba- 
sis of exact wavefunctions in the Coulomb field, as was 
done in Sec. IV.B.l.a, we now describe the solution ob- 
tained by Terekhov et al., 2008. These authors bypass 
the solution of the Dirac equation, obtaining instead an 
exact integral expression for the Green's function in a 
Coulomb field, using a proper-time approach common in 



QED (Mil'shtein and Strakhovenko, 1982). The main re- 
sult is that 



Sp(r) ^ -Q6{r) + Spdist 



(4.19) 



where Spdist{i') represents a positive charge distributed 
at r = oo (needed to satisfy the constraint of total zero 
induced charge). It is significant that this approach af- 
fords an exact expression for the dependence of Q upon 
g = Za, which is shown in panel (a) of Fig. 21. A series 
expansion of this dependence yields the following: 



Q(5)« 2 5 + 0.783 5^ + 1.398 5' 



(4.20) 



with each term corresponding to successive orders in per- 
turbation theory. The linear term is the one that ap- 
peared already in (4.17), at the RPA level. The next 
term in the expansion was also calculated perturba- 
tively by Biswas et ai, 2007. Interestingly, even though 
this problem is analogous to conventional QED vacuum 
polarization of a point charge, the perturbative coeffi- 
cients in Q{g) are not small, and increase with order, in 
stark opposition with the behavior known in 3D QED 
(Brown et ai, 1975). This offers another perspective 
upon the uniqueness of electron-electron interactions in 
graphene, for, even though the problem is on the surface 
analogous to the QED situation, the physics can be qual- 
itatively different. In this particular case, the difference 
seems to arise from the 2D dimensionality of the problem 
and the absence of Lorentz invariance in graphene, which 
renders the Coulomb interactions instantaneous. 

Inspection of the curve Q{g) in Fig. 21(a) reveals that 
it reaches 1 at g = 0.49, slightly before gc- This im- 
plies that, for a monovalent impurity (Z = 1), the non- 
interacting result predicts complete shielding before gc, 
insofar as Zcs{Z, a) = Z — Q{g) — > 0. Such strong renor- 
malization of the potential source immediately begs the 
consideration of interaction and correlation effects. They 
can be incorporated at the Hartree level by solving the 
self-consistent equation 



Zcfta = Za — aQ{Zcsa), 



(4.21) 



which encodes an infinite summation of a selected set of 
bubble diagrams (Terekhov et al., 2008). Since Q{g) is 
obtained exactly, one obtains the renormalized effective 
potential strength, Z^sa, with an accuracy much beyond 
the RPA. In addition, the reduction of Z^s with respect 
to the bare Z means that gc is also self-consistently renor- 
malized to (jc — Z^^a. The effect is shown in Fig. 21(b), 
which reveals that, as gc > gc, self-consistent screening 
delays the supercritical threshold because the condition 
ZcffCK = 0.5 requires a higher bare Z . This phenomenon 
is most striking for Z = 1 , in which case the supercritical 
point disappears altogether [gc < 1/2 even as Z — > oo), 
whereas g^='^ = 1.136 and g^='^ = 0.798. The predic- 
tion of this self-consistent Hartree renormalization of Z^tt 
would then be that impurities with Z = \ can never be- 
come supercritical. In addition, Hartree screening is suf- 
ficient to suppress the tendency for over-shielding of the 
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Coulomb center: as seen in the inset of Fig. 21(b), Z^ff 
remains always positive. 

An alternative approach to the Hartrec screening con- 
sists in treating the induced charge in linear response, 
Sp{q) = ^^(9)n(q), but taking into account electron- 
electron interactions perturbatively, via the renormal- 
ization of the coupling constant (Biswas et al., 2007). 
This is valid for small a (weak interaction), and leads 
to a result formally equivalent to (4.19), but where 
<^Pdist now arises from the electronic correlations. The 
distributed charge in the interacting case also has an 
anti-screening sign, but decays as l/r^, while the non- 
interacting (Spdist('') is zero everywhere, except at infin- 
ity. 

Even though the above considerations pertain to un- 
doped graphene, since all screening charge accumulates 
completely within a narrow distance, finite densities are 
not expected to alter the picture for as long as = Aakp 
remains large compared to the lattice scale a. 



2. Strong Coupling {g > Qc) 

In Sec. IV.B.l.c Hartree screening was shown to renor- 
malize gc and delay the critical threshold. Two important 
questions naturally arise: (i) since the self-consistent so- 
lution of (4.21) is uncontrolled, how certain can one be 
that the critical regime is reachable at all? (ii) So far we 
looked only at screening from the undercritical side (i.e. 
as long as Zcsa < 1/2). How can one address screening 
from the supercritical side, given that this regime cannot 
be reached perturbatively? 

The answer to these questions is far from trivial. In 
QED it is related to the ground state and stability of 
super-heavy nuclei {Z > 170), when the bound spectrum 
dives into the positron continuum (Fig. 18). Despite hav- 
ing received considerable attention throughout the 1970- 
80's (Greiner et ai, 1985), the fact that these systems 
require such high Z's, has turned it largely into an aca- 
demic problem. The exciting prospect about graphene is 
that impurities with Z = 1, 2 might already display su- 
percritical physics, in which case it would afford a bench- 
top test of some yet untested QED predictions. 

The essence of the difficulties in treating the super- 
critical regime clearly lies in its non-perturbative nature. 
Graphene, being gapless, is even more pathological be- 
cause of the infinite quasi-spectrum that appears in the 
hole channel [Fig. 18] . This quasi-spectrum is akin to an 
atom filled with infinitely many electrons and, as known 
from studies of heavy atoms (Landau and Lifshitz, 1981), 
it requires full consideration of correlations and interac- 
tions, and self-consistent techniques like Thomas-Fermi 
(Fermi, 1927; Thomas, 1927). 



a. Non-interacting Induced Charge. In Sec. IV. A. 3 we 
saw some unusual consequences for the DOS and cross- 
sections extracted from the exact solution of the Dirac 



equation for g > Qc- Now we address the correspond- 
ing induced charge obtained using the same procedure 
as in Sec. IV.B.l.a. Consideration of the exact wave- 
functions (Pereira et al., 2007) or the exact phase-shifts 
(Shytov et al., 2007a) leads to the conclusion that the 
supercritical partial waves contribute with an induced 
charge oc 1/r^. This could be expected on dimensional 
grounds: 5{r) and 1/r^ are the only dimensionally con- 
sistent possibilities in the absence of any intrinsic length 
scale in massless graphene. The exact induced density 
per partial wave reads (Shytov et al., 2007a) 

5n,{r) = ^V5^, (4-22) 

and, like the undercritical contributions, has a screening 
sign. The full induced charge is obtained from Sp{r) = 
-\e\5n{r), n(r) = J2\3\<gJ^j+T,\]\>gJnj' and has the 
general form 

6n{r) = SgA^ + Bsg6{r). (4.23) 

If 1/2 < g < 3/2 eq. (4.23) reduces to 6n{r) = 
{■Kg/2)5{r) + 2sgyjg'^ — gl/ {■n'^r'^). The general behavior 
(4.23) is also confirmed numerically by exact diagonaliza- 
tion of the tight-binding Hamiltonian in the honeycomb 
lattice, whose results are plotted in Fig. 20(b). 

b. Supercritical Protection. Unlike the undercritical 
regime, the additional power law decay in (4.23) causes 
a modification of Coulomb's law at large distances. But 
since we have a quasi-atom with all levels (4.8) filled, 
the non-interacting result in Eq. (4.23) cannot be the fi- 
nal answer. Each level is quasi-localized on the lattice 
scale, and should contribute significantly to shield the 
Coulomb center. For g not too much above gc we can 
follow an argument advanced by (Shytov et al., 2007a) 
that assumes electrons at some distance r feel the effect 
of a point charge consisting of the impurity subtracted 
from all the accumulated screening charge up to r. In 
other words, we introduce a distance dependent impu- 
rity strength, Zes{r) = Z — J^Sn{r)dr, and substitute 
(4.23) for Sn{r): 

^e.(r) = Z-f,-V^l,g^ (4.24) 

Since the log term represents the renormalization com- 
ing from screening at distances away from the center, we 
should replace {g = Za) — >■ {Z^ga = gcs). This leads to 
a self-consistent renormalization of the coupling that can 
be written in an appealing RG fashion as dgcff/d\og{r) = 
"^a-y/g^--^. In this way, it can be immediately seen 
that the coupling ^eff will "flow" to the constant value gc 
within a finite distance [see also (Gupta and Sen, 2009) 
for a related renormalization procedure]. As such, irre- 
spective of the bare Z, the system self-consistently re- 
arranges itself so that electrons at large distances never 
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feel a supercritical effective coupling. The undercritical 
(stable) situation is therefore protected. This reasoning 
agrees with expectations for the corresponding problem 
in QED, where it was shown that, within Thomas-Fermi, 
the vacuum polarization charge in super-heavy nuclei be- 
haves in such a way as to reduce Z to the threshold value 
(MiiUer and Rafelski, 1975). 

This is quite different from a metal, to the extent that 
graphene always leaves an universal amount of charge 
(Zc = gdoL) unscreened at large distances. Such behav- 
ior derives from the sharp transition between the under 
and supercritical regimes. On the one hand, the system 
wishes to screen as much charge as it possibly can. But, 
on the other, it cannot screen if 5 < ^c, therein lying the 
compromise that makes screening stop when Z reaches 
Z,. 



precisely, it applies ior XjZ <^ a <^ XjsfZ^ and pro- 
vided that log(r/i?) < l/a. Otherwise, for intermedi- 
ate electron-electron coupling [a ^ 1), the asymptotic 
screened potential should follow Fcft ~ ■^ce^/(eo''), with 
'Z'c = QcjoL = l/(2a). This result embodies the under- 
critical protection discussed above in Sec. IV.B.2.b, in- 
sofar as the supercritical core is always self-consistently 
screened so that Zoff — > Z^- Moreover, within the su- 
percritical core region, r < 2Za^R, the effective poten- 
tial decays as cx l/r'^/^. This obtains treating graphene 
as an ideal classical metal, under the assumption of 
quasi-complete screening in the core region (Fogler et at, 
2007). 



3. Finite IVIass 



c. Nonlinear Thomas- Fermi and Beyond. While the above 
approach is valid in principle only for g ^ gc, the fact that 
qualitatively supercritical graphene resembles a super- 
heavy atom suggests the use of TP theory, which is 
exact for atoms with Z ^ 00 (Lieb, 1981), and af- 
fords an approximation from the opposite limit g ^ gc- 
If wc wish to calculate how Coulomb's law is modi- 
fied in this regime we can calculate the total potential 
Feff(r) = V{r) + SV{r), where SV{r) = g / j^^dr' is 
the potential induced by the screening charge. Within 
TP we replace (5n(r') = n[fj, — V{r)] — n{ii), and the ho- 
mogeneous density depends on /j, via n = SE^J^ / {t^v^)- 
Solution of the resulting integral equation leads to the 
correction to Coulomb's law, which asymptotically reads 
(Katsnelson, 2006) 



Veff(r) 



K=ff(r) 



eor \_l + 2Za'^\og{r/R) 
Z 



eQr{qsry 



1 - 2ZoP- \og{qsR) 



(4.25a) 
(4.25b) 



valid toy fj, = Q, r ^ R and /i 7^ 0, rqs ^ 1 respectively, 
where Qs = Aa^/v is the screening length (4.18). One 
notes that the overall space dependence is formally the 
same as the one obtained within RPA, both at zero and 
finite density. Hence the bracketed coefficients in (4.25) 
can be interpreted as a renormalization of the valence. 
The important difference is that, in the limit Z — > 00 of 
interest in the context of TP, the nominal valence Z dis- 
appears from Vcf[{r), which thus becomes universal (and 
undercritical). Hence, even for strong impurities one can 
formally use perturbativc expressions for the screened po- 
tential, corrected for this renormalization of Z . 

It is important to emphasize that, since at this stage we 
are concerned with screening and corrections to the in- 
duced charge coming from electron-electron interactions, 
g = Za is no longer the relevant parameter alone, but 
both Z and a (that controls the interaction) indepen- 
dently. Por this reason, Pogler et al., 2007 have argued 
that the result (4.25) is valid only for small a. More 



We now briefly address the differences expected in 
the screening properties of charged impurities in massive 
graphene. We shall consider only the undoped situation, 
and assume fi = —to, such that none of the bound levels 
(4.6) are occupied. 



a. WeaJf Coupling (g < gc). It is clear that at weak 
coupling one can directly rely on perturbativc results 
(Sec. IV.B.l), and obtain the induced density from 
5n{q) = — ZVo(q)n(q). H(^)(<7) has been calculated 
in (2.14), and simple substitution yields the following 
asymptotics: 



f5(r) 



(5n(r) 



Za { -X^ 



-X 



log- 

-3 



r ~ a — > 
a < r < Ac , 
r > Ac 



(4.26) 



where Ac = l/{m,v) is the Compton wavelength, and 
a the lattice parameter of graphene. The short distance 
term is the same as found in the massless case (4.14,4.17), 
which makes sense given that when r <^ Ac the sys- 
tem does not "feel" the mass yet. It has a screening 
sign. However as the distance increases screening is in- 
creasingly suppressed, first weakly up to Ac, and then 
strongly, beyond Ac- In fact, since here 5n{q = 0) = 0, 
we have exactly / 5n{r)dr = 0. The meaning of this 
is simple: the total induced charge is zero. The system 
cannot screen beyond r > Ac because it is essentially 
an insulator (or a semiconductor with ji in the middle of 
the gap) . Notwithstanding, unlike a conventional insula- 
tor, gapped graphene shows a novel screening behavior at 
short distances, reflected in the live dependence of Sn{r) 
on the distance up to Ac- 



b. Strong Coupling (g > gc). In gapped graphene, screen- 
ing in the supercritical regime is qualitatively easier to 
understand, at least when g ^ gc- If the first level 
has just merged inside to hole continuum, its effective 
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probability density, |^'c(r)p, remains exponentially lo- 
calized, as described in IV. A. 2. a. Invoking complete- 
ness of the set of single-particle states, one can eas- 
ily show that the non-interacting induced charge follows 
(Pereira e< ai, 2008a) 

fc(r)« |^e(r-)p+(5npoi(r), (4.27) 

where Sup^iir) « J2E<-m \XE{r)\'^ - Ix^ir)] represents 
the vacuum polarization (i.e.: the induced charge coming 
from the full set of plane wave states), and is the same 
quantity that obtains in RPA (4.26). Clearly, the contri- 
bution from the supercritical state alone makes Sn{r) in 
(4.27) highly localized within the Compton wavelength, 
Ac- For all purposes, this state screens like a bound 
state would, and consequently one expects the impurity 
valence to be reduced by one unity times the degeneracy, 
N, of the level. But since iV = 4, this would imply, for 
the experimentally significant cases of Z ^ \, a, tendency 
to over-screen the Coulomb center. This bring us again 
to the role of interactions. The above would be true in 
the limit of weak interaction a <C 1. But, in that case, 
the supercritical regime would require Z ^ 1, which is 
not feasible. In the end, if supercritical systems arc to be 
produced, electron-electron interactions should be strong 
which, besides requiring the computation of the vacuum 
polarization in strong-coupling, brings the question of the 
renormalization of the bound levels themselves (Lamb 
shift). This situation, however, is completely analogous 
to the problem of super-heavy nuclei in QED, and an ex- 
tensive account of its particular features and difficulties 
can be found in Greiner et ai, 1985. 

C. From Single to Many Particle Interactions 

Coupling to an external Coulomb field can be seen as 
the zero-th order approach to the full many body elec- 
tron interactions in graphcne. The decisive difference 
that leaves graphene apart from standard electronic sys- 
tems is the existence of the supercritical region, which, 
for the Coulomb field, has the peculiarities discussed so 
far. Since the coupling constant in vacuum is a « 2, one 
can justifiably ask whether supercritical effects carry to 
electrons interacting among themselves. After all, even if 
a simplification, from a reference frame moving with an 
electron the problem becomes an impurity one again. 

1. Interacting Two Body problem 

The two particle problem has traditionally provided 
valuable insights into the full many-body phenomena 
in condensed matter [e.g. the Cooper pairing (Cooper, 
1956)]. The chiral nature of the electronic states, how- 
ever, precludes the usual decoupling between center-of- 
mass and relative coordinates, except for s-states in a 
quiescent center-of-mass (Sabio et ai, 2010b). Even so, 
these authors show that the supercritical collapse is a 



general effect present in the two body problem. In this 
case the critical coupling occurs at Uc = 1 and Uc = 2.24 
for s and p channels, respectively. The interacting two- 
body problem usually encodes much of the physics that 
the many-body system displays. One example is the 
study of pairing, pair condensation, and other processes 
which are dominated by two particle channel events. This 
has a clear relation with the issue of spontaneous gap 
generation, discussed in Sec. III.B. The prospect of exact 
solution of the two particle problem would afford more 
controllable means to explore this instability in graphcne. 



2. Excitons and Spontaneous Mass Generation 

It is noteworthy that the value ac = 1 quoted 
above is tantalizingly close to recent calculations of 
the critical coupling which precipitates a spontaneous 
mass generation and metal-insulator transition in un- 
doped graphene. Those values range from Uc ~ 0.8 
(Vafek and Case, 2008), to ac = 1.1 obtained within 
Monte Carlo (Drut and Lahde, 2009b) or by using the 
Schwinger-Dyson equation (Khveshchenko, 2009). As de- 
scribed at length in Sec. III.B, this metal-insulator tran- 
sition in graphene has been ascribed to the emergence of 
an excitonic instability beyond Uc- 

Recently the excitonic problem has been considered 
vis-a-vis the supercritical instability of the Coulomb 
center. Instabilities in the particle-hole channel ap- 
pear at critical couplings consistent with the above 
(Gamayun et ai, 2009; Wang et ai, 2010). For exam- 
ple, Gamayun et ai, 2009 show that solving the Bethe- 
Salpeter equation in graphcne leads to instability-prone 
tachyonic states (£'^ < 0) at ac = 1.6. Such states are 
the analogue in the two channel many-body language of 
the quasi-bound resonances for supercritical impurities, 
and a glimmer of supercritical effects in the fully inter- 
acting problem. 



D. Supercritical Physics in Experiments 

The non-perturbative nature of supercritical Coulomb 
impurities, and the associated analytical difficulties, pre- 
clude unequivocal predictions regarding the possibility of 
crossing the supercritical threshold. Experimental inves- 
tigation of this problem requires the ability to vary the 
strength of the Coulomb impurity and/or the electron- 
electron interactions. Control over the dielectric environ- 
ment provides a handle to tune interactions and impurity 
strength at the same time, via selection of eq- Experi- 
ments in this vein have been performed by Jang et at, 
2008 and Ponomarenko et ai, 2009, showing that it is 
possible to controllably tune the value of e by explor- 
ing substrates with different dielectric properties. Vari- 
ation of Z is a more delicate issue. Chen et ai, 2008 
have devised a way to add monovalent ions to graphene 
via K-irradiation, in quantities that can be controlled 
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with some precision. But exploration of the supercritical 
regime might require higher valences. For real impurities, 
the valence is determined by the nature of the impurity 
atom and the host system, and cannot be changed. One 
can, in principle, use ions of different valence, but here 
the difficulty lies in the fact that valences higher than 
Z = 2 are very unlikely. One possible alternative to this 
constraint imposed by nature, would be to resort to sharp 
STM tips, whose strong local field could mimic a strong 
local charge. As mentioned in the beginning of this chap- 
ter, the experimental exploration/confirmation of the su- 
percritical state would be rather important a milestone. 
Not only in understanding the physics of graphene, but 
because it would afford a glimpse to what might happen 
in the more fundamental QED situation. 

V. STRONG CORRELATIONS IN GRAPHENE 
A. Mass gaps in the honeycomb lattice 

Graphene is a semi-metal (SM) with gapless quasipar- 
ticles. The Dirac points in graphene are protected by 
the combination of sublattice and translational symme- 
tries of the honeycomb lattice. The point group symme- 
try of the honeycomb lattice, Cgi;, can be decomposed 
into the point group of the triangular sublattice and the 
Z2 sublattice symmetry group, C31, (g) Z2. Violation of 
sublattice symmetry leads to the opening of a mass gap 
in the Dirac Hamiltonian. This broken symmetry can 
be physically implemented cither by the Semenoff gap 
(Semenoff, 1984), which is induced by a staggered scalar 
potential that breaks the sublattice inversion symmetry, 
as previously discussed in Eq. (3.40), or by the Haldane 
gap (Haldane, 1988), where there is an additional broken 
time reversal symmetry (TRS) induced by the inclusion 
of circulating current loops with zero magnetic flux per 
unit cell, corresponding to a staggered magnetic field. In 
particular, a system that breaks inversion and TRS is 
susceptible to a "parity" anomaly, where the application 
of an electric field generates a net axial current flowing 
between the two valleys in graphene (Jackiw, 1984). 

In the presence of mirror symmetry along the z- 
axis, the spin-orbit interaction in graphene has the form 
(Kane and Mele, 2005) 

-Hso = Aso 51*1^^0® ^3® sa^-k,^, (5.1) 

k,(T 

where A50 is the spin orbit coupling gap, and S3 is 
the diagonal Pauli matrix in spin space. The other 
matrices follow the convention in the Dirac Hamil- 
tonian (2.6). The spin-orbit interaction in graphene 
breaks the spin degeneracy in the valleys, giving rise 
to spin polarized currents that flow along the edge 
states of the system — a quantum spin Hall state 
(Kane and Mele, 2005). Although the spin orbit cou- 
pling gap in graphene is rather small, A50 ~ 10~^meV, 
(Huertas-Herno et ai, 2006; Min et ai, 2006; Yao et ai, 



2007), it can be drastically enhanced either by curva- 
ture effects (Huertas-Herno et ai, 2006), or by impu- 
rities (Castro Neto and Guinea, 2009). The spin-orbit 
coupling is also logarithmically enhanced by Coulomb 
interactions (Kane and Mele, 2005), as discussed in Sec. 
III.B. When the mirror symmetry is broken either by a 
substrate or external electric field, an additional Rashba 
term is allowed 

'Hr ^ Xr *L^r3 ® (cTi ®S2-<72® si)*k,<T , (5.2) 

k,(T 

where A/^ > is the Rashba coupling. The induced gap 
is 2{Aso — Xr) for A_r < Aso, closing to zero when 
Xr> Xso (Kane and Mele, 2005). 

Kekule lattice distortions (Hou et ai, 2007), which 
break the translational symmetry of the lattice, also lead 
to the opening of gaps in graphene, whereas lowering 
the rotational symmetry of the Csv group, by stretching 
the honeycomb lattice in one direction, does not. In the 
presence of topological defects in the order parameter, 
such as vortices, the midgap states which are bounded 
to them allow the emergence of excitations with frac- 
tional statistics under vortex exchange (Chamon et at, 
2008a,b; Hou et ai, 2007; Seradjeh and Franz, 2008). In 
the superconducting case, the vortex core may sustain 
a quantum Hall state in the presence of a strong Zee- 
man coupling of the electrons with the magnetic field, 
which lifts the spin degeneracy (Herbut, 2010). In the 
most general case, where any spin, valley and pairing 
symmetries are allowed, 36 different types of instabilities 
that generate mass gaps in graphene have been classified 
(Ryu et ai, 2009). 

B. Charge and magnetic instabilities 

Although no evidence of mass gaps has been found in 
graphene, numerical results have predicted a semi-mctal- 
insulator (SM-I) transition in the presence of strong cor- 
relations. Quantum Monte Carlo (QMC) calculations on 
the Hubbard model for the honeycomb lattice at half 
filling predicted the opening of a Mott gap above the 
critical ratio U/t> 5 (Martelo et ai, 1997; Paiva et al., 
2005; SoreUa and Tosatti, 1992), where t « 2.8 eV is 
the hopping energy and U is the on-site electronic re- 
pulsion. A more recent QMC calculation has found 
a gapped AF state at half filling for U/t > 4.3, pre- 
ceded by an intermediate coupling insulating phase for 
3.5 < U/t < 4.3, which has been attributed to a gapped 
spin liquid state formed by short-range resonating va- 
lence bonds (Meng et al., 2010). An insulating antifcr- 
romagnctic (AF) ground state has been also predicted 
above U/t > 4 (Furukawa, 2001; Martelo et al, 1997). 
Variational (Hanish et ai, 1995) and mean field calcula- 
tions (Peres et al. , 2004) predicted the possibility of Na- 
gaoka ferromagnetism (where the polarization is maxi- 
mal) above a critical coupling both in the half filled and 
in the doped regimes. Although the validity of the Hub- 
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bard model in graphenc may be questioned sinee it does 
not include long range Coulomb interactions, it could be 
in principle justified if one accounts for strong screening 
effect from a substrate which can deplete the long range 
part of the interactions (or also, perhaps, by account- 
ing for dynamical screening effects from graphene itself 
(Reed et ai, 2010)), leaving only the short-range part of 
the electron-electron interactions. The extent of validity 
of the Hubbard model in graphenc is a subject of ongoing 
debate. 

The bare spin polarization in graphene is a 2 x 2 tensor 
(Peres et ai, 2004), 



n+-(q,r) = (5+(q,T)5-(-q,0)), 



(5.3) 



where 5+ and are the spin raising and lowering oper- 
ators in the two sublattices, x ~ a,b. Written in terms of 
the Green's function (2.11) with additional spin labels, 

k,s,s' — ± 



(5.4) 



where = 1+sk-cr/k, andi?s,CT(k) — sw|k|—// describes 
the two branches of the spectrum near the Dirac points. 
Since Ila^a = Hh b and ITa^f, = ^ by the honeycomb 
lattice symmetry, the eigenvalues of the spin polariza- 
tion are Hf/af ~ n+~ ± l^^tl ' '^hich correspond to 
ferromagnetic {+) and AF (— ) states. In RPA, the spin 
susceptibility is x = [1 ^ U&^^]~^&''-\ and the critical 
Hubbard coupling required for a divergence in the spin 
susceptibility in graphene is (Peres et ai, 2004) 



(5.5) 



The ferromagnetic transition translates in the condition 
= 2/p{(x) K, D"^ which is the Stoner criterion, 
where p{E) is the DOS and D the band width. The AF 
transition occurs at U^^ ~ D^/{D — |/i|). 

The application of an in plane magnetic field, B, splits 
the spin degeneracy at the Dirac points, creating two 
Fermi surface (FS) pockets with opposite spins. In- 
cluding the Zeeman coupling, Hb = X^o- "'^'^k.cr into 
the Hamiltonian, the spin polarized energy spectrum is 
Es.a{k) = sv\k\ + aB — /i. The nesting between the 
two Fermi surface sheets can produce a logarithmic di- 
vergence in the spin polarization in the limit \B\ ^ 
max(T, \fj,\) (Bercx et al, 2009), 



hW(o)^p(B) In 



\B\ 



max(r, 1^1) 



(5.6) 



This instability brings the possibility of a canted AF state 
in graphene. In the presence of Landau level quantiza- 
tion due to the application of an out of plane magnetic 
field, electronic interactions may lead to the formation 




FIG. 22 Semi-metal (SM) insulator transition predicted by 
the renormalization group analysis of the extended Hubbard 
model, in large A*' expansion. U is the on-site Hubbard cou- 
pling and V is the nearest neighbor site repulsion. Uc stands 
for the critical coupling. AF: antiferromagnetic phase; CDW: 
charge density wave state (Herbut, 2006). 



of quantum Hall ferromagnetic states at integer values of 
the filling factor (Nomura and MacDonald, 2006). The 
magnetic field has been also proposed as a source of a 
charge density wave (CDW) Peicrls distortion in the zero 
Landau level in graphene, breaking the parity symmetry 
between the valleys (Fuchs and Lederer, 2007). For a 
discussion of interaction effects at strong magnetic fields, 
see sec. VIII. 

For Dirac fermions in 2-1-1 dimensions, a CDW insta- 
bility translates into the phenomenon of chiral symmetry 
breaking (CSB), with spontaneous generation of a mass 
term that breaks the sublattice symmetry. The AF state 
is favored by strong on site repulsion and competes with 
the long range part of the Coulomb field, which can fa- 
vor either strong coupling fcrromagnctism (Peres et at, 
2004) or else excitonic CDW instabilities at strong cou- 
pling (Drut and Lahde, 2009a,b; Khveshchenko, 2001a;b; 
Khveshchenko and Leal, 2004; Liu et al, 2009). 

At large N, with N the number of fermionic flavors, 
the continuum limit of the Hubbard model in the honey- 
comb lattice falls in the universality class of the Gross- 
Neveu model (Gross and Neveu, 1974) for massless Dirac 
fermions in 2-1-1 dimensions, with four-fermion contact 
interactions. The extended version of this model accom- 
modates the short range piece of the Coulomb interaction 
involving the repulsion between nearest neighbor sites, V 
(Herbut, 2006). In addition to the Gaussian fixed point, 
which controls the semi-metal (SM) phase, the RG flow of 
the extended model was shown to be controlled by two 
other fixed points at large TV: an AF fixed point, and 
a CDW fixed point, both unstable towards the Gaussian 
fixed point at weak coupling, and having a runaway direc- 
tion to strong coupling when U or V are sufficiently large. 
The two fixed points compete, resulting in the phase di- 
agram shown in Fig. 22. The fact that the AF fixed 
point has only one unstable direction to leading order in 
1/A^ motivated the conjecture that the SM-I transition to 
the AF state is continuous and of the Gross Neveu type 
(Herbut, 2006). The symmetry analysis of the possible 
quartic terms has been discussed by Herbut et ai, 2009. 

The 1/N results were confirmed qualitatively by nu- 
merical renormalization group (NRG) calculations for 
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the extended Hubbard model in the honeycomb lattice 
(Raghu et at, 2008). In the presence of next- nearest 
neighbors repulsion, the NRG calculations suggested the 
possibility of competition between the CDW and spin 
density wave (SDW) phases with non-trivial topological 
insulating states, such as the quantum spin Hall (QSH) 
state, where TRS is spontaneously broken (Raghu et al., 
2008). Functional renormalization group (FRG) calcula- 
tions for the t — J model on the honeycomb lattice with 
on site and nearest neighbor repulsion also suggested the 
possibility of strong coupling CDW and SDW instabili- 
ties in graphene at half filling (Honerkamp, 2008). In the 
doped regime, the t — J model can favor the formation 
of superconducting states for J > 2t, either in the triplet 
or in the d-wave singlet channels (Honerkamp, 2008). 

In the high doping regime, the proximity of the Fermi 
level to the Van-Hove singularities, where the graphene 
DOS diverges logarithmically, may favor a Pomeranchuck 
instability (PI), rather than a gapped state. In that case, 
the redistribution of the electronic density generates a de- 
formation of the Fermi surface, which lowers the lattice 
Csv point group, instead of breaking the Z2 sublattice 
symmetry. In the extended Hubbard model at high dop- 
ing, the PI is favored by the repulsion between nearest 
neighbor sites, which renormalizes the kinetic energy at 
the mean field level, and competes with the on-site repul- 
sion, which favors a ferromagnetic state when the Stoner 
criterion is satisfied (Valenzuela and Vozmediano, 2008). 

When coated with metallic atoms that have a strong 
tendency to hybridize with the carbon orbitals, 
graphene can induce strong itinerant ferromagnetism in 
the metallic bands (Uchoa et ai, 2008b). 

C. Local magnetic moments 

For masslcss Dirac particles, the formation of local- 
ized states is usually harder than in usual Fermi sys- 
tems due to the Klein paradox, in which the fermions 
can easily tunnel through a barrier regardless of its 
height. Defects such as vacancies, where a carbon 
atom is knocked out from the plane, have been shown 
to generate localized states in graphene (Pereira et at, 
2006; Vozmediano et ai, 2005), and were recently ob- 
served in STM experiments (Ugeda et ai, 2010). Vacan- 
cies have also been found to host local magnetic states 
(Chen et ai, 2011; Yazycv and Helm, 2007). 

Short range interacting impurities can generate local 
resonances, which are quasi-localized states. At half- 
filling, the energy of the resonance, £o: is given by 
(Skrypnyk and Loktev, 2006; Wehling et ai, 2007) 

^° - So\n\sl/{D^--el)\ ' (''^^ 

where Uq is the scattering potential of the impurity and 
D is the bandwidth. The resonance induces accumula- 
tion of LDOS at the Fermi level around the impurity, 
p(r,uj), which decays as l/r (Bena and Kivclson, 2006), 



whereas the Friedel oscillations decay as 1/r^ for intra- 
cone scattering and as 1 /r for intercone scattering (Bena, 
2008). 

Besides defects, zigzag edges also lead to local mag- 
netism in the presence of interactions (for a more detailed 
discussion, see Sec. VI). In bulk graphene, a simple way 
to generate localized magnetic states is provided by the 
adsorption of adatoms with inner shell electrons. On 
the lattice, the adatoms can stay in different locations 
relative to the two sublatticcs in graphene. Transition 
metals are usually more stable sitting in the hollow site, 
at the center of the honeycomb hexagon (Chan et ai, 
2008), whereas simple molecules and atoms such as hy- 
drogen (H) tend to hybridize more strongly with the car- 
bons, sitting on top of them and generating a large local 
moment (Yazycv and Helm, 2007). In particular, H ad- 
sorption creates a midgap state (Boukhvalov et ai, 2008; 
Wehling et at, 2010c) and distorts locally the sp^ carbon 
bonds, which acquire sp^ character (Elias et ai, 2009). 
This distortion can induce a strong local enhancement 
of the spin-orbit coupling up to « 7 meV, as in dia- 
mond, and generate a strong local magnetic anisotropy 
(Castro Neto and Guinea, 2009). Adatoms can also 
form local moments from substitutional defects on single 
and double vacancies in graphene (Krasheninnikov et al., 
2009; Venezuela et ai, 2009). 

The heuristic criterion that describes the formation of 
a local magnetic moment is addressed at the mean field 
level by the Anderson impurity model (Anderson, 1961). 
In the top carbon case, assuming that the adatom sits on 
a carbon (see Fig. 23), say on sublattice B, the hybridiza- 
tion Hamiltonian is Hy = V ^„[flbcr{0) + h.c.], where fa 
(/J) annihilates (creates) an electron with spin a =t, i at 
the impurity. In momentum space, this translates into: 

nv = vJ2ifX.^ + blaL)- (5.8) 

p,cr 

If TT-o- — if a fa) is the occupation of the localized level for 




FIG. 23 (color on line) (a) Honeycomb lattice with an impu- 
rity atom. Black: sublattice A; White: sublattice B. Inter- 
section of the Dirac cone spectrum, E(k) = ±i)|k|, with the 
localized level Ef = £0: (b) eo > 0, (c) £0 < 0. 
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a given spin, the effeetive Hamiltonian of the level is 



(5.9) 



with £(j = Eq + Un^a, after a proper mean field decom- 
position of the Hubbard term Hu = Uf^f^fJ^fi, which 
accounts for the charging energy U to doubly occupy 
the level. The hybridized level becomes magnetic when 
^ n^. The occupation is derived self-consistently by 
integrating the /-electron DOS from the bottom of the 
graphene band up to the Fermi level /i, 



-Im 



(5.10) 



where T,ff{u}) is the self-energy of the localized electrons. 
In the cone approximation of the spectrum in graphene, 
for the top carbon case, 

^ff{io) =io[l^ Z-\io)] - iA\Lo\9{D - , (5.11) 

where A = ttV'^/D'^ is the dimensionless hybridization, 
D is the effective band width, and 



(5.12) 



gives the quasiparticlc residue, Z{uj), which vanishes log- 
arithmically at the Dirac points (w — > 0). 

Because of the vanishing DOS, the level 
broadening [given by ImI]^y(cL))] scales lin- 
early with the energy around the Dirac 
points (Gonzalez-Buxton and Ingersent, 1998; 
Skrypnyk and Loktev, 2006; Uchoa a/., 2008a; 
Zhang et ai, 2001). The DOS induced around the bare 
level, Ea, does not decay like a Lorentzian as in usual 
metals, but shows a long tail proportional to l/u. This 
tail induces several peculiar features in the magnetic 
states. For instance, a local moment is allowed to 
exist when the bare level is empty (eg < fi) or doubly 
occupied {eq + U > ji) (see Fig. 24). The presence of 
the Dirac point also breaks the symmetry around the 
line fJ. — Eq = U/2, and makes the scaling of the curves 
shown in Fig. 24 non-universal. Furthermore, there is a 
physical asymmetry between the cases where the level 
is above (eq > 0) or below {eq < 0) the Dirac point. 
When Eq = 0, as in the case of a vacancy, the level 
decouples from the bath and becomes magnetic for any 
/i > 0, regardless of the value of U (Pereira et ai, 2006; 
Uchoa et ai, 2008a). 

Since the chemical potential in graphene can be tuned, 
the formation of local magnetic states can be controlled 
by the application of a gate voltage (Uchoa et ai, 2008a). 
The low density of states around the localized level also 
makes the formation of local moments in graphene much 
easier than in usual metallic hosts. As a result the 
adatoms can achieve high magnetic moments at relatively 
small U (Cornagha et ai, 2009; Uchoa et ai, 2008a). 




5 10 
AD/U 



non-magnetic 



5 10 
AD/U 



FIG. 24 Boundary between magnetic and non-magnetic im- 
purity states in the scaling variables x = AD/U and y = 
(/i - eo)/U for eo > (a) and eo < (b). \eo\/D = 
0.029, 0.043, 0.029 and V/D = 0.14, 0.14, 0.04 for circles, 
squares and triangles, respectively. The upturn close to j/ = 1 
and a- — >■ on panel b) signals a crossover to the Fermi liquid 
regime /i, [/ ^ |eo| > 0, where the Dirac points are physically 
irrelevant. This feature is not visible in this scale when V is 
very small (triangles) (Uchoa et ai, 2008a). 



The formation of local moments is also affected by the 
specific location of the adatom in the lattice (Fig. 25). 
For instance, when the adatom sits in the center of 
the honeycomb hexagon (If-site), the tight-binding hy- 
bridization Hamiltonian is (Uchoa et ai, 2009) 



[Va,A{a,) + VbMi-^^)] /^(0)+h.c. , (5.13) 



where (i = 1,2,3) are the three nearest neighbor vec- 
tors of the honeycomb lattice, and V^^i {x = a, b) is the 
hybridization strength of the adatom with each of the 
nearest surrounding carbon atoms. In momentum repre- 
sentation. 



E 

per 



Vb^pbU fa + h.c. 



where 



1=1 



(5.14) 



(5.15) 



The top carbon case is recovered by setting T^.p = V and 
14, p = or vice-versa. For s-wave orbitals, Vx,i = V , 
whereas for in-plane /-wave orbitals the hybridization is 
anti-symmetric in the two sublattices, Va^i = —Vb,i = V. 
In the case of substitutional impurities (5'-sites), either 
Va,i = or Vb^i ~ 0. The quantum interference be- 
tween the different hybridization paths of the electrons 
can modify the energy scaling of the level broadening in 
Eq. (5.11) (Uchoa et ai, 2009), and can also change the 
shape of the Fano resonances in scanning tunneling spec- 
troscopy (STS) measurements, allowing a clear identifica- 
tion of the adatom position with an STS tip (Saha et al., 
2010; Uchoa et ai, 2009; Wehling et ai, 2010b). 
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D. Kondo effect 

The formation of a Kondo screening cloud around a 
magnetic moment is described by the Anderson Hamil- 
tonian (5.8) in the strong eouphng Hmit, [/ — >■ (X), where 
the valence fluctuations arc suppressed and the local mo- 
ment becomes a good quantum number. In the stan- 
dard mean field approach, the spin 1/2 fcrmionic fields 
are replaced by fermionic fields with larger degeneracy, 
N > m, which corresponds to an SU(iV) extension of 
the problem, with a corresponding Kondo Hamiltonian 
(Coqblin and Schrieffer, 1969) 

-Hk =JkJ2Y. ^Imflfm^^'.ra- , (5.16) 
mm' kk' 

where Jk ^ ^^/ko — mI is the Kondo coupling, -0^ 
are annihilation (creation) operators of the itin- 
erant electrons, and the local / fields arc constrained 
to a fixed occupancy. At the mean field level, which is 
asymptotically exact at large iV, the Kondo order pa- 
rameter can be extracted either from the standard slave 
boson approach to the Anderson model (Coleman, 1983; 
Newns and Read, 1987), or else by an equivalent path 
integral approach starting from the Kondo Hamiltonian 
(5.16) (Read and Newns, 1983). 

The application of these methods to semi-metals with 
a vanishing DOS, p{uj) = pq\uj\''\ with r > 0, resulted in 
the prediction of a Kondo quantum critical point (QCP) 
at half- filling (/i ~ 0). In that case, a Kondo cloud is 
expected for Jk > Jk ~ ^/{PoD^)i below the Kondo 
temperature (Withoff and Fradkin, 1990) 

Tk ^ \Jk - JkI" , (5-17) 

where ly = 1/r, and D is the ultraviolet cut-off. Since the 
scaling dimension of the hybridization V in the Anderson 
model is dim[y] = (1 — r)/2, the case r = 1 acts as an 
upper critical scaling dimension in the problem, where 
the scaling is marginal (Vojta and Fritz, 2004). In the 
marginal case, the Kondo temperature may have an addi- 
tional logarithmic scaling with the coupling, upon imple- 
mentation of an ultraviolet cut-off smoothly connected to 
the metallic case (r = 0) (Cassanello and Fradkin, 1996). 
Away from half-filling, there is a crossover to the usual 




FIG. 25 Two adatom configurations in graphene: a) the 
adatom (red circle) sits on top of a carbon atom, and b) 
the adatom (blue circle) sits at the center of the honeycomb 
hexagon, hybridizing equally with the two sublattices. Red 
arrows: nearest neighbor vectors. 
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FIG. 26 Schematic phase diagram around the Kondo QCP 
at half filling: temperature vs Kondo coupling. LM: local 
moment phase, where the Kondo cloud is suppressed. In 
the critical LM phase, quantum critical fluctuations dominate 
(Ingersent and Si, 2002). 



Fermi liquid case in the weak coupling regime, Jk < Jki 
where (Withoff and Fradkin, 1990) 

Tk « pcMr'\D/iiY{l - J-k/Jk) - • (5.18) 

Further studies based on NRG techniques 
(Fritz and Vojta, 2004; Gonzalez-Buxton and Ingersent, 
1998; Vojta, 2001) predicted a variety of fixed points. 
At half-filling, in the particle-hole symmetric case, 
Eq = —U 12, the Kondo problem has a metallic Kondo 
screened fixed point at r = 0, which evolves into a strong 
coupling fixed point for < r < 1/2. In this case, the 
strong {Jk > Jk) weak coupling {Jk < Jk) regimes 
are separated by a symmetric quantum critical point 
(SCP), whereas for r > 1/2 the local moment remains 
unscreened for all initial values of the Kondo coupling 
(Chen and Jayaprakash, 1995). In the particle-hole 
asymmetric case {^ = 0, U ^ —2eo), for r > r* w 0.375, 
the weak and strong coupling regimes are separated 
by an asymmetric critical point (ACP). For r < r* , 
the particle-hole symmetry is dynamically restored 
(Fritz and Vojta, 2004; Gonzalez-Buxton and Ingersent, 
1998). 

The phase diagram around the QCP is schemati- 
cally shown in Fig. 26. The critical local moment fluc- 
tuations were studied by Ingersent and Si, 2002, who 
found linear uj/T scaling of the dynamical spin sus- 
ceptibility at the critical point for < r < 1. In 
the marginal case, r ~ I, there are logarithmic cor- 
rections to scaling (Cassanello and Fradkin, 1997). The 
Kondo problem for gapless excitations was also exten- 
sively studied in the context of magnetic impurities 
in d-wave superconductors (Borkowski and Hirschfeld, 
1992; Cassanello and Fradkin, 1996, 1997; Polkovnikov, 
2002; Polkovnikov et ai, 2001; Vojta and Bulla, 2001; 
Zhang et ai, 2001; Zhu and Ting, 2000). For a review, 
see Balatsky et ai, 2006. 

In the graphene case, where r = 1, the Dirac 
fermions in the bath have an additional pseu- 
dospin structure, which motivated several proposals for 
multichannel Kondo physics (Cassanello and Fradkin, 
1996; DeU'Anna, 2010; Sengupta and Baskaran, 2008; 
Zhu et ai, 2010). The Kondo resonance in graphene 
has been calculated with NRG by Cornaglia et al. , 2009. 



At half filling, the local DOS around the impurity can 
be spontaneously enhanced by the formation of midgap 
states due to the scattering potential of the impurity 
(Hentschel and Guinea, 2007), frustrating the Kondo 
QCP. 

At finite doping, the Kondo temperature has an 
exponential dependence with the DOS at weak cou- 
pling, allowing the Kondo cloud to be tuned by gat- 
ing (Sengupta and Baskaran, 2008). In the crossover 
regime, at J = Jc, the scaling of the Kondo tempera- 
ture with doping becomes power law, Tk oc Recent 
NRG calculations in graphene have found a particle-hole 
asymmetric scaling of the Kondo temperature with dop- 
ing, Tk oc where a; = 1 for /i > and x = 2.6 
for /X < (Vojta et al., 2010), in contradiction with 
the mean field and poor man scaling analysis for the 
marginal case (Vojta et ai, 2010). In the presence of 
Landau level quantization, the Kondo temperature has 
reentrant behavior as a function of the chemical poten- 
tial (Dora and Thalmeier, 2007). 

Looking at the problem on the lattice, ab initio calcu- 
lations on Cobalt have found that the interplay of spin 
and orbital degrees of freedom can give rise to an SU(4) 
Kondo effect in graphene when the spin orbit coupling is 
strong enough (Wehling et ai, 2010a). Another ab initio 
calculation accounting for dynamic correlations, also on 
Co, has identified the possibility of a spin 3/2 Kondo 
effect, involving multiple orbitals (Jacob and Kotliar, 
2010). From a tight-binding perspective, for a spin 1/2 
impurity, the hybridization Hamiltonian (5.14) can be 
written in the diagonal basis 



(5.19) 



-± p.cr 



where c±,ko- = (l/%/2)[6ko- ± ((/>k/l</'k|)akcr] are the 
fermionic operators that diagonalize the graphene Hamil- 
tonian (2.2), 0k = s'*^' is the tight-binding hop- 
ping matrix clement defined by Eq. (2.3), and a = ± 
labels the conduction and valence bands. O is a phase 
factor, which accounts for the symmetry and position 
of the localized orbital with respect to the sublattices 
(Uchoa et ai, 2009), 



V2V 



aV* — ? 

m 



where Vx^p is the hybridization as defined in Eq. (5.15). 

As in metals, the Anderson Hamiltonian in graphene 
can be mapped into the spin exchange Hamiltonian by 
a canonical transformation (Schrieffer and Wolff, 1966). 
In the large U limit, the spin exchange Hamiltonian be- 
tween the magnetic adatom and the graphene electrons 
is (Uchoa et al, 2011) 



kk' oca' 



cL',^',k''^Cc<,a,k, (5.21) 




0.4 0.6 

FIG. 27 Kondo coupling vs. chemical potential in graphene 
for U = 1 eV and V = 1 eV. The Kondo coupling can be 
controlled by gate voltage across the weak (J <C Jc) and 
strong coupling (J > Jc) Kondo regimes, where Jc is the 
critical coupling at half-filling. 



where cr = (cri , (T2 , ) are the Pauli matrices, S 
^flcrf^i is the localized spin, and 



j(m) 



(eo - M)(fo + U - ^i) 



<0, 



(5.22) 



is the exchange coupling defined at the Fermi level, 
/i. Within the tight-binding description, we realize 
that the determinant of the exchange coupling matrix 
in Eq. (5.21) is identically zero, det[JaQ'] = 0, and 
hence the exchange Hamiltonian (5.21) can be rotated 
into a new basis where one of the hybridization chan- 
nels is decoupled from the bath (Pustilnik and Glazman, 
2001). The eigenvalues in the new diagonal basis are 
JuMM' = JJ2a^*a,k^aM' ^ud Jy = 0, implying that 
the one-level exchange Hamiltonian (5.21) maps into the 
problem of a single channel Kondo Hamitonian, Jig = 
■ Sk,k'j where s is the itinerant spin, in 
spite of the implicit valley degeneracy. A multi-channel 
description of the one-level problem is nevertheless possi- 
ble for example in graphene quantum dots, in the contin- 
uum limit, where valley and angular momentum channels 
become good quantum numbers. 

Unlike the situation in metals, the exchange cou- 
pling in graphene can be controlled by gating 
(Jacob and Kotliar, 2010; Uchoa et ai, 2011), as shown 
in Fig. 27, in particular when the chemical potential is 
brought to the proximity of the localized level, where 
the Kondo coupling becomes resonant. This effect opens 
the possibility of tuning J to the vicinity of the criti- 
cal coupling that sets the crossover between the weak 
and strong coupling regimes. In this region, at finite 
doping, quantum criticality is reminiscent of the frus- 
trated QCP at n = 0. Since the width of the Kondo 
peak in the spectral function is set by the Kondo tem- 
perature only, the gating effect permits measuring the 
quantum critical scaling of the Kondo temperature with 
doping (Uchoa et ai, 2011; Vojta et ai, 2010) directly 
with STM probes (Saha et ai, 2010; Uchoa et ai, 2009; 
Wehling et ai, 2010b; Zhuang et ai, 2009). 
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E. RKKY interaction 

The Ruderman-Kittel-Kasuya-Yosida (RKKY) inter- 
action between two local spins is obtained by integrat- 
ing out the itinerant fermions in Eq. (5.21), which gives 
Hrkky = —J^Xij Si ■ where Xij is a two point corre- 
lation function, with i, j indexing the positions of the lo- 
cal spins. In momentum space (Brey et al., 2007; Saremi, 
2007; Uehoa et al, 2011), 

v-noi = V M^y /K'(k + q)]-/[g.(k)] 



where (omitting the aa' labels for simplicity) 



k.q 



p)*x p.y p,x p.*v 

^a,k^Q,k^Q',k+q^Q',k+q - 



(5.23) 



(5.24) 



with x^y = A, B, H, S etc, indexing the position of the 
spins on the lattice, i?Q,(k) = a|0k| — Mi f is the 
Fermi distribution. A^^q = -^kq = 1/4 for spins on 
the same sublattice whereas 



1 , 0k<?!>k 
—(y.ct 

4 |0k||0k+q| 



^k-Hq 



(5.25) 



for spins on opposite sublattices. In the continuum limit, 
where the spectrum is linearized around the Dirac points. 



M 



AB 
k,q 



jaa'e^^^-^+'i, where 6 is the angle between k and 
k + q (Brey et a/., 2007). 

At half-filling, kp = 0, the Fermi surface collapses into 
points and the RKKY interaction is mediated by inter- 
band transitions, which polarize the vacuum as in QED. 
In this case, the Friedel oscillations disappear and the 
sign of the interaction is ferromagnetic for spins on the 
same sublattice and anti-ferromagnetic for spins in op- 
posite sublattices (Brey et al., 2007; Saremi, 2007). In 
the overdoped regime, at = i, the nesting among 
the Van Hove singularities in graphene reverses the sign 
of the RKKY interaction compared to the fi = case 
(Uchoa et ai, 2011). 

At long distances, the spatial decay of the 
RKKY is when ^ is at the neutrality point 

(Brey et ai, 2007; Cheianov and Fal'ko, 2006; Saremi, 
2007; Vozmediano et ai, 2005; Wunsch et al, 2007). 
Away from half filling, the Friedel oscillations are re- 
stored by the intraband transitions and the RKKY in- 
teraction decays at r ^ I/Zcf as l/r^, similarly to the 
2DEG case (Brey et al, 2007; Wunsch et ai, 2007). For 
H or S site spins formed in C'^y symmetric orbitals, the 
RKKY interaction decays with a fast power law 1 /r^ at 
half filling (Uchoa et ai, 2011). In carbon nanotubes, 
the RKKY interaction decays as 1/r for top carbon 
spins and as for H site spins in isotropic orbitals 
(Kirwan et ai, 2008). 

When distributed regularly on top of graphene, mag- 
netic adatoms such as hydrogen (H) can form macro- 
scopic magnetic states at room temperature (Zhou et ai, 
2009). In the disordered case, H atoms in particular can 
cluster on top of graphene due to rippling. On top of 



a ripple, the sp'^ carbon (C) bonds are spontaneously 
stretched by the curvature and acquire sp'^ character. 
Contrary to the perfectly flat case, the adsorption of H 
atoms on top of the hills helps to stabilize the ripples 
(Boukhvalov and Katsnelson, 2009). The interplay be- 
tween the correlations due to the ripples and the RKKY 
interaction among the H spins can generate magnetore- 
sistance hysteresis loops and a variety of magnetic spin 
textures (Rappoport et ai, 2009). 



F. Superconductivity 

The observation of proximity induced supercon- 
ductivity in graphene junctions has stirred a lot 
excitement in the field of mesoscopics (Heersche et al, 

2007) . The Dirac nature of the quasiparticlcs gives 
rise to ballistic transport on a micron scale and allows 
graphene to sustain supercurrents in long junctions, 
the size of the coherence length in the superconduct- 
ing metallic leads (Du et al., 2008; Heersche et ai, 
2007; Miao et al., 2007; Ojeda-Aristizabal et al., 2009). 
The experimental realization of the proximity ef- 
fect motivated theoretical studies of the differential 
conductance (DC) in normal-superconductor (NS) 
interfaces in graphene (Beenakker, 2006; Burset et ai, 

2008) , graphene nanoribbons (Rainis et ai, 2009), and 
in graphene normal-insulator-superconductor (NIS) 
junctions (Bhattacharjee and Sengupta, 2006). Due 
to the Dirac nature of the spectrum, at half-filling, 
the Andreev conversion of an electron into a hole at 
the interface between a normal and a superconduct- 
ing region involves specular refiection rather than 
retro reflection (Beenakker. 2006). The specular 
Andreev reflection leads to the presence of Andreev 
modes in SNS junctions that propagate along the 
graphene edges at the interface with the superconductor 
(Titov et ai, 2007). The Josephson current in graphene 
SNS junctions was studied by Titov and Beenakker, 
2006, followed by Bergman and Hur, 2009; 
Maiti and Sengupta. 2007; Moghaddam and Zarcyan, 
2006, and Black- Schaffer and Doniach, 2008. Pos- 
sible applications involving the proximity effect 
in graphene include proposals for valley sensors 
(Akhmerov and Beenakker, 2007), current switches 
(Linder et ai, 2008; Lutchyn et ai, 2008), and a spin 
current filter (Greenbaum et ai, 2007). A review on 
Andreev and Klein tunneling processes in graphene can 
be found in Beenakker, 2008. 

These experimental developments in transport moti- 
vated a surge of interest in the possibility of making 
graphene an intrinsic superconductor. Graphene par- 
ent compounds, such as the graphite intercalated ma- 
terials CaCe and KCg, are low temperature supercon- 
ductors, although neither graphite nor alkaline metals 
alone superconduct (Csanyi et ai, 2005; Hannay et at, 
1965; Weller et at, 2005). Even though intrinsic su- 
perconductivity has not been observed in the single 
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FIG. 28 Superconducting order parameter Aij = Aie'*^ 
{j = 1, 2, 3), with phases along tire tliree different bond di- 
rections in the lattice. 



layer so far, a few different superconducting mecha- 
nisms have been proposed. One possibility is a plasmon 
mediated mechanism in graphene coated with metallic 
adatoms, in which the plasmons of the metallic band 
mediate the attraction between the graphene electrons 
(Uclioa and Castro Neto, 2007). When isolated islands 
of metallic atoms arc adsorbed on top of graphene, su- 
perconductivity can also be induced by proximity ef- 
fect (Feigel'man et ai, 2008). Another possibility is the 
Kohn-Luttinger mechanism, which explores the proxim- 
ity of the Fermi surface to the Van-Hove singularities in 
the high doping regime (Kohn and Luttinger, 1965). In 
this scenario, the superconductivity can be mediated by 
a purely electronic mechanism, when the interactions be- 
come attractive along a specific direction of the BZ near 
the Van-Hove singularity (Gonzalez. 2008). The super- 
conductivity can also be mediated by in plane or out 
of plane flexural phonons (Lozovik and Sokolik, 2010). 
In graphene, strong doping regimes can be currently 
achieved by chemical adsorption of alkaline metals, such 
as potassium (Griineis et ai, 2009; McChesney et al., 
2010; Uchoa et ai, 2008b), or with metal contacts 
(Giovannetti et al, 2008). 

Alternative proposals include edge state supercon- 
ductivity, induced by the large DOS at the edges 
(Sasakia et at, 2007), or strong correlations, which so far 
have not been observed in graphene. As in the cupratcs, 
the antiferromagnetic attraction between spin singlets 
on nearest neighbor sites has been proposed as a pos- 
sible pairing channel in graphene, provided the on site 
Hubbard repulsion is strong enough to suppress the lo- 
cal fluctuations (Pathak et ai, 2010). Gonzalez et ai, 
2001 considered the possible competition between ferro- 
magnetic and superconducting states in graphene sheets 
through a renormalization group analysis accounting for 
Coulomb interactions. A recent functional renormaliza- 
tion group calculation has proposed the possibility of a 
strongly correlated SDW state that gives way to a singlet 
superconducting instability in the c?-wave channel, or else 
a CDW solution that allows a triplet pairing instability 
in the /-wave channel (Honerkamp, 2008). In two-layer 
graphene, the possibility of excitonic pairing of electrons 
in one layer with holes in the other one has been consid- 
ered (Kharitonov and Efctov, 2008; Min et al, 2008). 



Regardless of the microscopic origin, the superconduct- 
ing state in graphene can be analyzed based on the sym- 
metries of the order parameter in the honeycomb lattice. 
On the lattice, the electrons in graphene carry spin, angu- 
lar momentum and sublattice quantum numbers. There 
are four possible pairing channels: singlet/triplet spin 
channels, and same/opposite sublattices. In the singlet 
case, if we restrict the analysis to nearest neighbor site 
interactions only, two competing order parameters can 
be identified: 

Ao = goiar^aji) = go{hi^bji), (5.26) 

which corresponds to an s-wave state, and Ai, defined as 

^lAj = gi{a.i^bji - aiibj^) (5-27) 

for nearest neighbors and zero otherwise, where go and 
gi are the coupling strengths. In momentum space, the 
latter state is described by 

3 

Ai,k = ^Ai,,e''^'-'^, (5.28) 

i=l 

where Ai^i = Ai(aj) are the real space pairing ampli- 
tudes along the three different bond directions in the 
honeycomb lattice (see Fig. 28). In the simplest case 
the pairing amplitudes are the same, Ai_i = Ai, and Ai 
is real, giving 

Ai,k = Ai0k, (5.29) 

where = X]i=i e'*' '^' gives the hopping matrix el- 
ement in the single particle tight-binding spectrum 
(Uchoa and Castro Neto, 2007). This order parameter 
represents the pairing between electronic states in oppo- 
site sides of the BZ, and preserves all the physical sym- 
metries of the honeycomb lattice, including point group 
and time-reversal symmetry, Ai k = A|[ j,, where the 
momentum k is measured with respect to the center of 
the BZ, at the F point. In real space, this order param- 
eter (OP) has extended s-wave symmetry. If expanded 
around the Fermi surface centered at the Dirac point K, 
from the perspective of the quasiparticle excitations near 
the Fermi level, 

Ai,K+p = Aic*^ (p, + ipy) (5.30) 

describes a p + ip state in one valley and p — ip in the 
opposite one (Uchoa and Castro Neto, 2007). This state 
is therefore a p + ip state with additional valley degener- 
acy. Unlike the case of conventional p + ip superconduc- 
tivity, the time reversal operation involves an additional 
exchange of valleys, preserving the TRS of this state, and 
we shall refer to it as p -I- ip. 

Another possible paring symmetry is the state 
(Black-Schaffer and Doniach, 2007; Jiang et al., 2008) 

Aij = Aie^t^'^/^b" ^ (5 31) 
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FIG. 29 (color online) Order parameter (OP) amplitude, 
|Ai,k|, in the BZ: (left panel) Ai,j = Ai with j = 1,2,3 
indexing the three different bond directions of the crystal [see 
Eq. (5.28)] and (right panel) Ai^- = e'^""^^^, which describes 
a flux phase. Light colors represent higher amplitude. Dirac 
points are located at the K points, at the edges of the BZ. 
In all dark spots, the OP has p + ip symmetry around the 
respective high symmetry points. In the three light spots on 
the right panel, the OP has s-wave symmetry around the K' 
points. 



j = 1,2,3, which describes on the lattice a real space 
pairing wavefunction with d.j.2_y2 + id^j^-wave symme- 
try, breaking TRS. This broken symmetry is caused by 
the circulation of plaquette current loops, which amounts 
to global circulation of current along the edges. The 
low energy description of this state around the Dirac 
points is a combination of s-wave in one valley and 
p + ip state in the opposite valley (Jiang et ai, 2008), 
as shown in Fig. 29. At the mean field level, this state 
was shown to have lower energy than the pure p + ip 
state (Black-SchafFer and Doniach, 2007). Due to the 
broken TRS, disorder and quantum fluctuations, which 
are paramount in a 2D system, may strongly inhibit the 
coherence of the d + id state. Other alternatives are the 
degenerate states with d^^-y-i and dxy-^S''vc symmetries, 
represented by the Ai , pairing amplitudes (2,-1,-1) 
and (0, 1, —1), respectively (Black-SchafFer and Doniach, 
2007). These states conserve TRS but lower the crystal 
point group symmetry. 

In the spin triplet channel, the OP is a superposition of 
Sz = — 1,0,-t-l states. Since on-site pairing is forbidden 
by the Pauli principle, for nearest neighbors interaction 



the triplet superconducting states are A*^ 
with cr =t,| for Sz = ±1, and A; 



in the 5^ = channel. The OP in this case is commonly 
defined as a 2x2 tensor, 



(5.32) 



where the Pauli matrices act in spin space, and d.y = 
— dji is an anti-symmetric tensor, violating parity. The 
case where the OP d has a single vector compo- 
nent describes the spinless fermionic case, discussed by 
Bergman and Hur, 2009. The possibility of spin triplet 
states beyond nearest neighbors in the Sz = Q channel 
was recently examined in a variational cluster approxi- 
mation calculation (Sahebsara and Scncchal, 2009). An- 
other possibility is a Kekulc superconducting state in the 
triplet channel, which breaks the translational symmetry 
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FIG. 30 Phase diagram between the s-wave and effective p-\-ip 
phases in the spin singlet channel. On the left: ji = case, 
which is quantum critical. Right: ^ 7^ case. Continuous 
lines represent second order transitions, and dashed lines rep- 
resent flrst order transitions (Uchoa and Castro Neto, 2007). 



of the lattice and allows the presence of topological exci- 
tations (Roy and Herbut, 2010). 

At the level of nearest neighbor sites, the electron- 
electron interaction can be decomposed into an effective 
local Hubbard term, 



and a non-local part, 



^ b- 



(5.33) 



(5.34) 



In the singlet pairing channel, the non-local term can 
be decomposed into Hj = gi {~^ij^i3 + ^Ij^y) ' 

plus one body terms that can be absorbed into the 
chemical potential fi. Dij = Oit^ji ^ '^'i^it ^ stan- 
dard singlet pair operator and Bij = o,ia^j<^ ^ 
bond operator. Decomposition of the interaction at 
the mean field level with (Bij) ~ results in the 
graphcne tight-binding Hamiltonian for the supercon- 
ducting phase, H'^ ~ J^k ^'k'^^'I'k + Eq, where 



i?o = -|Ao|V5o-3A?/gi, 



and 



-tct>l 
AS 



\* 

1,-k 

A* 



Ao Ai,k \ 
Ai,_k Ao 



(5.35) 



(5.36) 



is the Bogoliubov-de Gennes matrix written in the sub- 
lattice and Nambu basis 5'k = (flkt, ^kt, f^-ki' ^-k^)- 

The Hamiltonian (5.36) can be diagonalized in a basis 
of Bogohubov quasiparticles: H'' = J2kas -^k,Q,s^k,a.s + 
Eq, where is the quasiparticle number operator and 
s,a = ±1. In the isotropic case, Ai_k = Ai(/)k, the spec- 
trum is E)f^^a,s = cx.E]^^s, with (Uchoa and Castro Neto, 
2007) 



(i|0k| + snY + (|Ao| + sAi|0k|) , (5.37) 
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FIG. 31 Dependence of the gap, normalized by the band cut- 
off D, on g in the weak (g < gc) and strong coupling {g > gc) 
sectors for ^ = (red line), ^i/D — 0.1, and 0.3. The model 
has a QCP at half filling (Uchoa et al, 2005). 

where the phase of the OP Aq is locked in with Ai, and 
Al is real. The electronic gap described by the spectrum 
(5.37) is 

Eg^2\tAo-fiAi\/^Jt^ + Al. (5.38) 

In the p + ip state (Aq = 0, Ai ^ 0), Eg is proportional 
to the deviation of the chemical potential away from half- 
filling, and a.t fi = this state becomes quantum criti- 
cal and gapless. The instability in this case translates 
into the renormalization of the Fermi velocity, where 
t = ty^l + Af is the renormalized hopping amplitude, 
instead of the opening of a gap (Uchoa and Castro Neto, 
2007). Minimization of the free energy 

F = -T ^ In [2 + 2cosh(Fk,s/T)] + Eq , (5.39) 

k,s 

with respect to Aq and Ai gives a set of two cou- 
pled BCS-like equations, and leads to the phase diagram 
shown in Fig. 30. At half-filling, = 0, the emergence 
of superconductivity is controlled by quantum critical 
lines in the parameters go and gi, with critical values 
gjj = — TTii^/D and g^ = —Attv^/D^, in the linear cone 
approximation, where D is an ultraviolet cut-off and v 
is the Fermi velocity near the Dirac point (Castro Neto, 
2001; Marino and Nunes, 2006; Uchoa and Castro Neto, 
2007; Zhao and Paramekanti, 2006). For finite n, there 
is a crossover to the standard Fermi liquid case at weak 
coupling, as shown in Fig. 30. 

When Ai,j = Aie'(27./3)j jggg (5.28)], the elec- 
tronic wavefunctions collect different phases along the 
different bond links, which gives rise to a current fiow, 
and the d + id state cannot coexist with an isotropic TRS 
s-wave state. The gap properties of the d + id state and 
the differential conductance in SN junctions were derived 
by Jiang et ai, 2008. The Josephson current for this 
state in SNS junctions was calculated by Linder et al, 
2009. 

In the s-wave state (we assume Aq to be real), the 
gap variation with the coupling at half filling, near the 
quantum critical point — — ttu^/Z?, is (Castro Neto, 
2001) 

Ao = D{1 - g^go) . (5.40) 



Away from half-filling, the gap crosses over to 
(Uchoa et ai, 2005) 

Ao - 2|Ai|exp [D{1 - g^o/9o)/\t^\ - 1] (5.41) 

for ^ Ao, which corresponds to the weak coupling 
BCS limit, where g g^, a.s shown in Fig. 31. The 
1^1 /Ao ^ 1 limit corresponds to the strong coupling 
regime (g > gc), and the intermediate coupling region 
near g ^ gc sets the crossover scale between the two 
regimes at finite fi. Non-equilibrium effects in the pres- 
ence of a dissipative environment may also lead to a dis- 
sipation driven quantum phase transition away from half 
filling (Takei and Kim, 2008). 

At mean field level, the critical temperature at ^ = is 
Tc = Ao/21n4, whereas in the opposite limit, \fi\ ^ Ao, 
Tc = 'jAo/tt, as in the BCS case, where In 7 « 0.577 
is the Euler constant (Uchoa et ai, 2005). Of course, in 
two dimensions there is no true long range order. The su- 
perconducting transition is of Kosterlitz-Thouless (KT) 
type and coherence is actually lost at much lower tem- 
peratures due to the role of thermal fluctuations, which 
unbind vortex and anti- vortex pairs above the KT tran- 
sition temperature, at Tkt < Tc- The mean field result 
indicates the onset of critical fiuctuations where the am- 
plitude of the Cooper pairs is completely destroyed, al- 
though the phase coherence is suppressed much earlier, 
at Tkt- The KT fiuctuations of the SC order parame- 
ter have been considered by Loktev and Turkowski, 2009, 
without accounting, nevertheless, for the chiral nature of 
the quasiparticles in graphene. 

Zero field thermodynamic properties, such as the spe- 
cific heat at fixed volume, Cy = —T{d^F/dT^)v, can be 
extracted from the free energy (5.39). For an isotropic 
condensate of Dirac fermions, the jump of the specific 
heat at the phase transition, normalized by the specific 
heat on the normal side, is (Uchoa et ai, 2005) 

SCv = 2(ln4)V[9C(3)] « 0.35, (5.42) 

at half-filling. In the |^|/Ao ^ 1 limit, the jump grows 
to the standard BCS value 6Cv ^ 12/[7C(3)i w 1.43. 

The Mcissncr effect in graphene, which describes 
the expulsion of an external magnetic field by the 
circulation of diamagnetic supercurrents, has been 
recently examined by Kopnin and Sonin, 2008 and 
Uchoa and Castro Neto, 2009. In the presence of vor- 
tices, the Bogoliubov de Gennes equations for Dirac 
fermions in 2-1-1 dimensions allow the presence of zero 
energy modes (Jackiw and Rossi, 1981) which are bound 
to the vortex cores. For a vortex with vorticity n (the 
winding number of the OP), Ao ~ |A„(r)|e*"''^, with 
(r, (j)) as cylindrical coordinates. The physical solutions 
allowed by the boundary conditions at the center of the 
vortex and at infinity result in n zero modes at half filling 
(Ghaemi and Wilczek, 2012). The subgap spectrum and 
the wavefunctions in the vortex core have been derived 
by Bergman and Hur. 2009; Seradjeh, 2008. Away from 
half filling, for odd vorticity n, there is only one energy 
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branch that crosses zero energy for zero angular momen- 
tum. For n even, no subgap branch intersects zero energy, 
and no exact zero modes exist (Bergman and Hur, 2009; 
Khaymovich et ai, 2009). Because of the fermionic de- 
generacy in the valleys, the topological zero modes do not 
lead to fractional statistics under vortex exchange, as in 
conventional p+ip superconductors, unless additional in- 
teractions that lift the fermionic degeneracy are included 
(Herbut, 2010). Vortex zero modes for excitonic conden- 
sates in bilaycrs have been discussed by Seradjeh et at, 
2008. 



VI. INTERACTIONS AT BOUNDARIES AND LATTICE 
DEFECTS 

A. Surface states 

The vanishing density of states of graphene at the neu- 
trality point implies that localized states can exist at the 
Dirac energy, much in the same way as localized states 
appear inside a forbidden energy gap in semiconductors 
and insulators. In order for these states to be normal- 
izable, special boundary conditions are required. These 
conditions imply the breaking of the translational sym- 
metry of the lattice, so that they can only exist near 
edges or defects. 

The most extensively studied examples are the sur- 
face states which exist at graphene zigzag edges, where 
the lattice is abruptly terminated (Fujita et ai, 1996; 
Nakada et ai, 1996). Such edges have been observed 
in graphene flakes (Girit et ai, 2009; Jia et ai, 2009), 
and also in graphite (Niimi et al., 2005). As the local- 
ized states form an energy band of zero width, the lo- 
cal density of states at the Dirac energy near a zigzag 
edge changes from zero to infinity, and the electron com- 
pressibility becomes divergent. Interactions of arbitrarily 
small strength lead to instabilities when the Fermi energy 
lies at the Dirac point. A mean field analysis showed that 
a short range Hubbard interaction can lead to a ferromag- 
netic ground state (Harigaya, 2001; Harigaya and Enoki, 
2002). In zigzag ribbons with two edges, the spins at the 
two edges are aligned antiferromagnetically, see Fig. 32. 
These early theoretical results, based on the tight bind- 
ing approximation, were later confirmed by calculations 
based on the Local Density Approximation (Pisani et at, 
2007; Son et ai, 2006). The ferromagnetic order re- 
mained when the dangling bonds at the zigzag edges 
where saturated by hydrogen, which probably is closer to 
the actual experimental situation. The optimization of 
the atomic positions at the edges leads to reconstructed 
phases with gaps, where the spin up and spin down bands 
do not overlap near the gap, suggesting a half metallic 
phase (Son et ai, 2006). Other phases, with ferroelec- 
tric properties (Fernandez-Rossier, 2008) or canted mo- 
ments have been studied (Jung and MacDonald, 2010). 
A sketch of the magnetization induced near a zigzag edge 
of a graphene ribbon is shown in Fig. 32. Recent experi- 




FIG. 32 (Color online). Sketch of the magnetization at the 
zigzag edges of a graphene ribbon. 



ments (Enoki and Takai, 2009; Joly et ai, 2010) confirm 
the existence of magnetic moments at graphene edges. 

The effects of the electron-electron interaction on the 
midgap states has also been studied beyond the mean 
field approximation. The calculations show that the fer- 
romagnetic phase is stable when the band of localized 
states is half filled. Both a local onsite interaction or 
the long range exchange effect lead to this phase. At 
very low fillings, electrons tend to form a charge density 
wave state, similar to a Wigner crystal (Wunsch et al., 
2008a,b). More complex correlated states are possible 
at other fillings. The fact that the midgap states at a 
zigzag edge resemble the wavefunctions of Landau lev- 
els, in that the momentum parallel to the edge and the 
spatial extension are coupled, leads to the intriguing 
possibility of states similar to the Laughlin wavefunc- 
tions which describe the Fractional Quantum Hall Effect 
(Wunsch et ai, 2008a). 

At long distances, straight graphene edges of arbitrary 
orientation other than armchair can support midgap 
states, as zigzag edges (Akhmerov and Beenakkcr, 2008). 
Hence, local magnetic moments can be a generic prop- 
erty of abrupt graphene edges. Zigzag edges and vacan- 
cies in bilayer (Bernal) graphene also give rise to midgap 
states, at least when only the direct nearest neighbor in- 
terlayer hopping is included (Castro et ai, 2008a), and 
magnetic moments can be formed at the edges of bi- 
layer graphene (Sahu et ai, 2008). Models which include 
other interlayer hoppings lead to sharp resonances near 
edges and vacancies. These results suggest that moder- 
ate interactions can produce local moments in graphene 
bilayers or in three dimensional graphite. The combi- 
nation of the Zeeman field associated with magnetic or- 
dering, and the spin orbit coupling can lead to phases 
characterized by quantized spin currents at the edges 
(Soriano and Fernandez-Rossier, 2010). 
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FIG. 33 (Color online). Sketch of the magnetization induced 
near a vacancy. 



B. States at vacancies and cracks 

Midgap states can occur in other situations where the 
translational symmetry of the lattice is broken. Similarly 
to the case of surface states at a zigzag edge, interactions 
will lead to the spin polarization of these states. The 
simplest situation where the existence of a partially lo- 
calized midgap state can be demonstrated is a lattice 
vacancy (Pereira et al., 2006, 2008b). This analysis can 
be extended to multilayer samples (Castro et al., 2010). 

The existence of these states has been confirmed by 
STM spectroscopy on vacancies in irradiated graphite 
(Ugeda et ai, 2010). It can be expected that interactions 
lead to the formation of a magnetic moment around the 
vacancy. The formation of local moments near vacan- 
cies is consistent with the observation of ferromagnetism 
in irradiated graphite (Barzola-Quiquia et ai, 2007; 
Chen et ai, 2011; Esquinazi et ai, 2003; Ohldag et ai, 
2007; Ramos et ai, 2010). Absorption of hydrogen leads 
to similar effects to those of a vacancy, including the 
formation of magnetic moments (Yazyev, 2008). Other 
dopants, like carbon atoms and NO2, also lead to the 
formation of spins (Lehtincn et at, 2003; Wehling et at, 
2008b). 

A sketch of the magnetization induced near a graphene 
vacancy is shown in Fig. 33. The moment associated 
with the localized level around the vacancy is coupled 
to the extended states, leading to the possibility of the 
Kondo effect. Some differences between usual mag- 
netic impurities and the situations described here can 
be expected: i) The vacancy or adatom modifies sig- 
nificantly the electronic density of states, rendering in- 
valid perturbative treatments which relate the magni- 
tude of the exchange coupling to the unperturbed elec- 
tronic structure. The phase shift induced in the conduc- 
tion band remains significant, even near the Dirac energy 
(Hentschel and Guinea, 2007). ii) The localized state is 
orthogonal to the extended states. Hence, the coupling 
between the local moment and the conduction band does 
not take place via virtual hops between the two types 
of states. Instead, it can be expected that the electron- 
electron interaction favors a ferromagnetic alignment of 
the local moment and the spins of the conduction elec- 



trons. 

Spins at different vacancies interact ferro- or antifer- 
romagnetically (Brey et ai, 2007; Palacios et ai, 2008), 
depending on whether the vacancies occupy the same or 
different sublattices. At half filling, the RKKY interac- 
tion mediated by the tt band decays as l/|r — r'p, and 
it goes to the l/|r — r'p dependence typical of a two 
dimensional electron gas at finite carrier concentrations 
(Cheianov and Fal'ko, 2006). Voids or cracks can be con- 
sidered an intermediate case between vacancies and edges 
(Vozmediano et ai, 2005). They also support localized 
spins at the boundaries. 



C. Midgap States and Random Gauge Fields 

Midgap states in bulk graphene can also be induced by 
magnetic fields (see below), or by strains which mimic 
the effect of a magnetic field (Guinea et ai, 2008b). 
These states have been analyzed using the tight bind- 
ing approximation (Guinea et ai, 2008b), or by means 
of the Local Density Functional method (Wehling et al., 
2008a). Corrugations and wrinkles also induce midgap 
states in graphene (Katsnelson and Prokhorova. 2008; 
Pereira et ai, 2010). The presence of these states en- 
hances the effects of the interactions. Mean field calcula- 
tions suggest the formation of magnetic moments, which 
will order ferro- or antiferromagnetically (Guinea et at, 
2008a,b). 

A random strain distribution leads to a random 
gauge field acting on the electrons. The changes 
in the electronic density of states induced by a ran- 
dom gauge field have been studied by RG techniques 
(Horowitz and Doussal, 2002; Ludwig et ai, 1994). Re- 
lated problems arise at the transition between plateaus 
in the Quantum Hall Effect, and in d-wave supercon- 
ductors. It can be shown that, above a certain disorder 
strength, a random gauge field leads to a divergent den- 
sity of states at the Dirac energy (Horowitz and Doussal, 
2002; Riu and Hatsugai, 2001). This divergence leads to 
a vanishing electron compressibility, and enhances the ef- 
fects of interactions in the same way as the midgap states 
considered earlier. A random gauge field, A(r), can be 
characterized by a dimensionless number. A, 

(A^(r)A,(r')) - A5^,(5(2)(r - r'). (6.1) 

If the gauge potential is assumed to arise from random 
corrugations of average height h and length then A ~ 
h'^/{a^P), where a is the lattice constant (Guinea et ai, 
2008a,b). A similar parameter can be defined if the gauge 
potential is due to topological defects, such as disloca- 
tions (Gonzalez et ai, 2001). The regime A 1 corre- 
sponds to ripples large enough to accommodate midgap 
states, leading to a divergence in the density of states. 
The changes in the density of states induced by a gauge 
field can be written as a logarithmic rcnormalization of 
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FIG. 34 (Color online). Sketch of the magnetization induced 
at the edges of a quantum dot. 



the Fermi velocity 



V ^ V 



1 - cA loe 



where c is a numerical constant, and A is a high mo- 
mentum cutoff of the order of the inverse of the lattice 
constant. 

The scaling towards lower Fermi velocities in eq. 6.2 
can be combined with the RG analysis of the long 
range Coulomb interaction (Foster and Aleiner, 2008; 
Foster and Ludwig, 2006a,b; Staubcr et ai, 2005). Dis- 
order tends to increase the density of states near the 
Dirac energy, while interactions lead to the opposite ef- 
fect. To lowest order, this analysis leads to a line of fixed 
points characterized by a finite disorder and finite inter- 
actions, as discussed in Sec. III.A.l, see Fig. 9. The tem- 
perature and frequency dependence of properties such 
as the conductivity or the specific heat acquire anoma- 
lous exponents (Herbut et ai, 2008). For high disorder, 
A > 1 , it can be shown that a gapped state is more stable 
than the gapless density of states expected in the absence 
of interaction effects (Guinea et ai, 2008a). 

Certain strain configurations lead to effects simi- 
lar to those induced by a constant magnetic field 
(Guinea et at, 2010). The possible ways in which the 
degeneracies of these states are lifted by the interactions 
have been studied (Herbut, 2008), and new phases, with 
properties similar to those of topological insulators may 
exist. It is worth noting that STM experiments sug- 
gest the existence of very large effective fields due to 
strains, -Be// ~ SOOT, in small graphene bubbles under 
high strains (Levy et ai, 2010). The effects of electron- 
electron interactions in this regime remain unexplored. 
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FIG. 35 (Color online). Single energy peaks and Coulomb di- 
amonds in a graphene quantum dot, see (Ponomarenko et ai, 
2008). 



VII. INTERACTION EFFECTS IN MESOSCOPIC 
SYSTEMS 

A. Magnetism in quantum dots 

Mesoscopic samples have a large ratio between the 
perimeter and the area. Midgap states localized at 
the edges can have a significant weight in the to- 
tal density of states, and interaction effects are en- 
hanced. Early calculations for planar carbon molecules 
(Stein and Brown, 1987; Tyutyulkov et ai, 1998) showed 
gaps associated with the electron-electron interaction, 
and magnetic moments at the edges. A large mag- 
netic moment can be found in triangular graphene flakes 
(Fernandcz-Rossier and Palacios. 2007), where the three 
boundaries have the zigzag orientation, and the carbon 
atoms at the edges belong to the same sublattice. 

As mentioned previously, edges of arbitrary orienta- 
tions, except the armchair direction, support midgap 
states (Akhmerov and Becnakker, 2008). Hence, local 
moments and magnetism can be expected in graphene 
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FIG. 36 (Color online). Sketch of the extension of edge states 
in a graphene quantum dot. 



quantum dots of any shape, provided that the termi- 
nation at the edges is abrupt. Model results suggest 
that this is the case, and the orientation of the mo- 
ments at the edges depends on the type of sublattice 
at the edge (Fernandez-Rossier and Palacios. 2007), as 
sketched in Fig. 34. Away from half filling, correlated 
states with unsaturated magnetization, and charge den- 
sity wave states are also possible (Romanovsky et ai, 
2009; Wunsch et ai, 2008a). The charging of a quantum 
dot leads to a substantial rearrangement of the electronic 
levels, in a similar way to the well studied orthogonality 
catastrophe in metals (Anderson, 1967; Wunsch et al., 
2008a). The conductance can acquire a non trivial volt- 
age or temperature dependence, as in a Luttinger liquid 
(Kane and Fisher, 1992). 

A simple estimate of the number of magnetic moments 
in a quantum dot can be obtained by assuming that the 
average density of edge states is of order pedge ~ c x 
[R/{aW)], where c 1 is a numerical constant, R is 
the radius of the dot, a is the lattice spacing, and W is 
the bandwidth of the band of edge states (Wimmer et al., 
2010). The Coulomb interaction within each state, which 
leads to the formation of local moments is Ec ~ e^/R x 
log(i?/a), see below. Naturally, one has to replace -^• 
e^/co in all formulas, but we do not write the dielectric 
constant explicitly in this section. The states which arc 
spin polarized are those whose distance from the Fermi 
energy is less than Ec- This condition, combined with the 
estimate for pedge , gives a maximum number of magnetic 
moments within the dot, N ss EcPedge ~ c x [e'^/{aW)] x 
log(i?/a). This number is not too large. For W ^ 0.3 — 
0.5eV and R - lOOnm we obtain iV - 10 - 20. The total 
magnetic moment of the dot depends on the sign of the 
couplings between the edge spins, see Fig. 34. 

Experimentally, there is evidence which suggests the 
formation of local moments in small graphene flakes, of 
dimensions 10 — 50nm (Sepioni et al, 2010). 



FIG. 37 (Color online). Sketch of a graphene ribbon with 
disordered edges as a series of quantum dots. 
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FIG. 38 (Color online). Graphene point contact coupled to a 
quantum dot, see (Stampfer et al, 2009). 



B. Charging effects. Coulomb blockade 

Graphene quantum dots of many shapes and di- 
mensions are being extensively studied (Avouris et at, 
2007; Bunch et ai, 2005; Gucttinger et al, 2008; 
Giittinger et ai, 2009; Han et ai, 2007; Huard et ai, 
2007; Molitor et ai, 2009a; Moser and Bachtold, 2009; 
Ozyilmaz et ai, 2007; Ponomarenko et ai, 2008; 
Stampfer et ai, 2008; Wilhams et ai, 2007). Sin- 
gle electron effects have been observed in many of 
them. Experiments show clear evidence of charg- 
ing effects in graphene quantum dots, as evidenced 
in the diamond patterns formed by the resonances 
in the conductance through the dot as a func- 
tion of gate and bias voltages (Guettinger et al., 
2008; Giittinger et ai, 2009; Molitor et al., 2009a,b; 
Moriyama et ai, 2009; Moser and Bachtold, 2009; 
Ponomarenko et ai, 2008; Ritter and Lyding, 2009; 
Schnez et ai, 2009; Stampfer et ai, 2008), see Fig. 35. 

The electrostatic interaction between electrons leads 
to Coulomb blockade, which modulates the energy dif- 
ference between levels, and induces non Ohmic features 
in the conductance through the dot. In a graphene quan- 
tum dot of dimension R, the electrostatic energy required 
to add a unit of charge scales as e^/R. The mean level 
spacing between extended states in a ballistic dot scales 
as v/R. As the dimensionless parameter a = e^/(eof) in 
graphene is of order unity, the energy scales associated 
with charging and confinement effects are comparable. 
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The edge states discussed earlier can lead to charging 
energies larger than those for extended states. Assuming 
that these states are delocalized along the perimeter of 
the ribbon, over a scale L ^ R, see Fig. 36, and width a 
comparable to the lattice spacing, see Fig. 36, the charg- 
ing energy becomes (e^/i?) x log(_R/a) (Wimmcr et ai, 
2010). 

Charging effects can also modify the transport prop- 
erties of narrow graphene ribbons. Irregularities in the 
edges may induce the formation of constrictions and 
quantum dots, as sketched in Fig. 37, where charging 
effects will lead to a transport gap. In a nanoribbon of 
width W, the typical size of these dots will also be W, 
and the transport gap will be of order e'^/W. In the 
absence of charging effects, a ribbon will have confined 
subbands, separated by gaps of order v/W. Hence, the 
similarity between the energy scales arising from quan- 
tum confinement and charging effects, which exists in a 
quantum dot, also exists in a graphene ribbon. An exper- 
imental realization of an all graphene circuit with a point 
contact coupled to a quantum dot (Stampfer et al., 2009) 
is shown in Fig. 38. This setup can be used to count the 
passage of charges through the quantum dot. 

Experiments in graphene nanoribbons are compatible 
with the relevance of charging effects (Han et al., 2010, 
2007; Todd et al., 2009). Some observations can be ex- 
plained by a model of dots formed in the ribbon con- 
nected through many channels with the rest of the struc- 
ture. Such a strongly coupled dot always shows Coulomb 
blockade effects, unless there is a perfect transmission 
through one or more of the channels. The effective 
charging energy, however, is strongly renormalized by 
the coupling between the dot and the rest of the system 
(Sols et al., 2007), ~ e'^/We'S, where g is the con- 
ductance, in dimensionless units, of the junction between 
the dot and the electrodes. In general, g ^ x kpW, 

where T is the transmission amplitude of a given channel. 

The electron-electron interactions can be studied in 
mesoscopic samples through their effect on the magneto- 
conductance at low magnetic fields. These experiments 
probe the phase coherence of electrons at low temper- 
atures. This quantum effect is suppressed due to the 
dephasing induced by the interactions. Electronic quan- 
tum coherence also gives rise to the universal conduc- 
tance fluctuations observed in disordered metals, which 
are also reduced by the dephasing due to interactions. 
The dephasing length shows a temperature dependence 
consistent with the expected behavior in a dirty metal, 
-^0 ~ (ff^^)/[2^1c)g(5)]: where g is the conductivity in 
dimensionless units (Tikhonenko et al., 2009) (sec also 
(Chen et al., 2010)). This dependence is replaced by a 
cx T^^ in high mobility samples (Tikhonenko et al., 
2009), as expected in a clean Fermi liquid. Experiments 
that tune the ratio between the dephasing length and 
the mean free path (Moser et al., 2010) show a variety of 
regimes, interpolating between weak and strong localiza- 
tion. 
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FIG. 39 (Color online). Sketch of the successive splittings of 
the Landau levels as the magnetic field is increased, a) Spin 
states are split first, and then the valley degeneracy is broken, 
b) Valley degeneracy is lifted first, followed by the breaking 
of spin degeneracy. 
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FIG. 40 (Color online). Splittings of the Landau levels in 
graphene as function of magnetic field, see (Zhang et al., 
2006). 



VIM. INTERACTIONS IN STRONG MAGNETIC FIELDS 

A comprehensive review of graphene in magnetic field 
has recently appeared (Goerbig, 2011), and here we only 
mention some of the main effects. The electronic en- 
ergy bands of graphene in a strong magnetic field col- 
lapse into Landau levels. In the absence of disorder. 
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FIG. 41 (Color online). Resistance of a suspended graphene 
sample as a function of carrier density for two different mag- 
netic fields. R. V. Gorbacliev, D. C. Elias, A. S. Mayorov, A. 
A. Zhukov, K. S. Novoselov, A. K. Geim (unpublished). 

the electronic compressibility diverges when the chemi- 
cal potential coincides with the energy of a Landau level, 
and the effects of the interactions are enhanced, as in 
other two dimensional metallic systems. The typical scale 
of the electronic wavcfunctions is the magnetic length, 
is = y/f^/{eB) = ^$o/(27rB), where B is the applied 
field and $o is the quantum unit of flux. The separation 
between levels is of order v/£b, while the relevant scale 
for interaction effects is /ts- 

There are two sets of Landau levels in graphene, one 
for each valley. In addition, graphene has the n = 
level, which combines electron and hole features. Hence, 
interactions can break either the valley degeneracy or the 
spin degeneracy. The long range part of the Coulomb in- 
teraction is independent of the valley index. The n = 
Landau level is localized in a given sublattice, and its 
degeneracy can be lifted by interactions which break the 
symmetry between sublattices, like the coupling to out of 
plane optical phonons (Fuchs and Lederer, 2007). Hence, 
the removal of the spin and valley degeneracies of the 
Landau levels due to interactions depends on other en- 
ergy scales (Goerbig, 2011), such as the Zeeman splitting, 
or the nearest neighbor repulsion, for the case n = 0. A 
sketch of the possible symmetry breaking patterns as a 
function of magnetic field is shown in Fig. 39. Early ob- 
servations of splittings between Landau levels are shown 
in Fig. 40, see (Zhang et ai, 2006). 

It is usually assumed that the Zeeman splitting 
is much smaller than the other energy scales. Cal- 
culations suggest that the spin degeneracy is lifted 
first, leading to excitations with combined spin and 
valley indices (Abanin et ai, 2007; Alicea and Fisher, 
2006; Goerbig et ai, 2006; Gusynin et ai, 2009; 
Nomura and MacDonald, 2006; Shibata and Nomura, 
2008; Wang et al, 2008; Yang et ai, 2006). The fourfold 



spin and valley degeneracy when the Zeeman coupling 
is neglected gives a new SU (4) symmetry, which may 
lead to new features, not observable in other two di- 
mensional electron gases (Goerbig and Regnault, 2007; 
Toke and Jain, 2007). The formation of Landau levels 
favors the excitonic transition which can also exist in 
the absence of a magnetic field (Gusynin et ai, 2006). 
The spin split n = level leads to spin polarized edge 
states (Abanin et al, 2006, 2007; Fertig and Brey, 2006; 
Shimshoni et al., 2009) where the orientation of the spin 
depends on the sign of the current, as in topological 
insulators (Hasan and Kane, 2010; Qi and Zhang, 2011). 

A magnetic field oriented parallel to the plane does not 
give rise to Landau levels. In neutral graphene, it leads 
to metallic states with electrons and holes polarized in 
opposite directions, providing another route towards an 
excitonic transition (Aleiner et ai, 2007). 

Experiments show that, indeed, the spin and val- 
ley degeneracies of Landau levels in graphene are 
lifted (Giesbers et ai, 2007, 2009; Jiang et al, 2007; 
Zhang et ai, 2006). The opening of a gap in the n = 
level in graphene has been extensively studied, and 
a metal insulator transition with critical features con- 
sistent with a Bcrezinskii-Kostcrlitz-Thouless transition 
has been reported (Amado et ai, 2009; Checkelsky et ai, 
2008, 2009). 

The most striking manifestation of the interac- 
tions in the presence of a strong magnetic field is 
the Fractional Quantum Hall Effect. Early theo- 
retical calculations showed that the FQHE could be 
stable in graphene (Apalkov and Chakraborty, 2006; 
Castro Neto et ai, 2006; Toke et al, 2006). The condi- 
tions for the FQHE are the existence of sharp Landau lev- 
els and sufficiently strong electron electron interactions. 
The analysis of FQHE states in graphene can be done 
in a similar way to that of a two dimensional electron 
gas. The main difference is a change in the pseudopoten- 
tials which describe the interactions between electrons 
in a given Landau level, because the wavcfunctions in 
graphene and in a two dimensional electron gas differ. 

This Fractional Quantum Hall Effect was extensively, 
but unsuccessfully, sought in samples deposited on Si02. 
Suspended samples, which showed a much higher elec- 
tron mobility, did not exhibit the FQHE, using the stan- 
dard experimental four terminal setup. The observa- 
tion of the IQHE in suspended bilayer graphene using 
a two terminal setup (Feldman et ai, 2009) led quickly 
to the discovery of the FQHE in single layer graphene 
(Bolotin et ai, 2009; Du et ai, 2009), using the same 
technique. More recently, four terminal measurements in 
high mobility suspended samples (Ghahari et ai, 2011), 
and also samples deposited on a new substrate, boron ni- 
tride (Dean et ai, 2011, 2010), also show the FQHE. In 
two terminal measurements, the existence of the FQHE 
is inferred from plateaus of the longitudinal resistance at 
carrier densities which correspond to fractional fillings of 
Landau levels, see Fig. 41. The v = 1/3 state turns out to 
be more robust than in other materials, like GaAs, which 
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exhibit the FQHE, and it can be observed at tempera- 
tures greater than lOK. Fractional plateaus dX v = 2/3 
and V ^ \/2 have also been reported. Theoretical cal- 
culations suggest that the so called Moore-Read ground 
state at fillings with even denominators, which leads to 
the existence of non Abelian anyonic quasiparticles, is 
not favored in graphenc (Wojs et al.. 2011). 



IX. INTERACTIONS IN BILAYERS 

Bilayers are the building blocks for 3D stacks of 
graphene, such as graphite. In a bilayer one has two 
parallel graphene sheets, separated by an equilibrium 
distance similar to the interlayer distance of graphite 
(3.35A) (Dressclhaus and Dresselhaus, 1981). The rel- 
ative position of the two graphene layers is not unique, 
and this leads to different stacking arrangements of the 
bilayer, and even more possibilities for multilayers, or 
graphite. The most stable configuration seems to be the 
so-called Bernal AB stacking, in which the two layers 
are rotated by 60°. As a consequence, one of the the 
sublattices in the lower layer (say, sublattice A) is ver- 
tically aligned with one of the sublattices of the upper 
layer (say, sublattice B) [see Fig. 42(a)]. Notice that 
this particular rotation leads to a breaking of sublattice 
symmetry between layers. As a first approximation, the 
electronic coupling between the layers can be described 
in terms of the hopping of electrons between the near- 
est neighbor atoms in different layers with an energy t±^ 
(also known as 71 « 0.39 eV in the graphite literature 
(Castro Neto et al., 2009a). Another possible arrange- 
ment between the layers is the fully aligned configuration, 
also called AA stacking. In both AB and AA stacking, 
the unit cell is comprised of 4 atoms, and has the same 
2D extension as the unit cell of a single layer; this im- 
plies that the Brillouin zone is precisely the same as in 
monolayer graphene. 

Notice, however, that these configurations are just a 
few of an infinite series of commensurate structures be- 
tween two layers, the so-called twisted bilayer graphene 
(Lopes dos Santos et al., 2007). The problem of com- 
mensurate and incommensurate structures always ap- 
pears when two crystalline materials are superimposed, 
as in the case of bilayers. For commensurate structures, 
the angle between the layers is not arbitrary but follows 
a well defined sequence (Lopes dos Santos et al., 2007). 
Obviously, different angles lead to different broken sym- 
metries and hence to different electronic states. When 
the angle of rotation is 60 degrees, as in the case of 
the Bernal structure, the sublattices arc noncquivalent, 
which leads to a broken sublattice symmetry and hence 
to a putative gap opening. For other angles, there is no 
broken sublattice symmetry but the unit cell is enlarged 
as the rotation angle becomes smaller. In this case the 
massless Dirac dispersion has to be preserved for sym- 
metry reasons (Li et al., 2010; Lopes dos Santos et at, 
2007; Mole, 2010). From this perspective, the Bernal con- 




FIG. 42 a) Top view of a graphene bilayer; white and black 
circles: top layer carbon atoms; gray and red: bottom layer, 
b) four-band spectrum of the bilayer, ±i5^(p), with 7 = ± as 
shown in Eq. (9.4), near the corner of the Brillouin zone, c) 
Brillouin zone with high symmetry points, d) Illustration of 
the four band spectrum around the K point. 

figuration is an exception. The twisted bilayer graphene 
presents a very rich physics of its own that we will not 
cover in this review. Instead, we will focus on the Bernal 
configuration which is the most studied case. 

We will start from the minimal tight-binding model 
for Bernal bilayers, which includes a basis with two ad- 
ditional layer fiavors (denoted by an ovcrbar), 

'I'k.cr = (ak,(T, &k,<T, bk.cr, flk.cr) , (9.1) 

with a i representing the spin. The resulting Bloch 
Hamiltonian is then a 4 x 4 matrix with two sublattice, 
and two layer degrees of freedom. 
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(9.2) 

where tj_ w 0.39eV is the interlayer hopping, and t w 
2.8eV is the in-plane, nearest neighbor, hopping ampli- 
tude. The momentum dependence is contained in (jjk, 
which is the same as for a monolayer (2.3). The band 
structure associated with eq. (9.2) consists of four non- 
degenerate bands given by 

Eik)=±^{t^± ^tl+U^\M'). (9.3) 

An expansion k = K + p around the K points of 
the BZ when ?;|p| <C t shows that the four-band tight- 
binding spectrum (9.3) resolves into four hyperbolic 
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bands (Nilsson et al., 2006), as shown in Fig. 42(b), and 
whose form reads: 



(9.4) 



with w w 6 eVA being the Fermi velocity (the same Fermi 
velocity of a monolayer), and 7 = ±1. The Bernal stack- 
ing explicitly breaks the sublattice symmetry in each 
layer, causing an energy split of t± between the two 
7 = ±1 branches, E+ and at p = (see Fig. 42). 
Due to a degeneracy at the K points, the two symmetric 
branches +E- and —E^ touch there, resulting in a gap- 
less spectrum. Just as in a monolayer, the Fermi surface 
of an undopcd bilaycr reduces to only two points, at K 
and K'] but now the valence and conduction bands have 
a finite curvature and, hence, notwithstanding the ab- 
sence of a gap, the effective electronic degrees of freedom 
are massive, but still chiral. The degeneracy at K is pro- 
tected by the Z2 symmetry between the two layers only 
(McCann, 2006) , and can be lifted with arbitrarily small 
perturbations, such as the ones induced by a bias volt- 
age, by polarizing the two sheets (Zhang et al, 2009), or 
else by independently changing the carrier concentration 
in each layer (Ohta et ai, 2006). This property opens 
the exciting prospect of using graphene bilayers as mate- 
rials with a gate-tunable band gap (Castro et ai, 2007; 
Castro Neto et ai, 2007; Min et ai, 2007). 

We stress that the low energy effective theory of bilay- 
ers remains Lorentz invariant, in the following sense. The 
rotation of 7r/3 between layers breaks the sublattice sym- 
metry leading to 2 pairs of massive Dirac particles at the 
K {K') point. Nevertheless, the system remains metallic 
because two bands, belonging to different pairs, touch in 
a point. More explicitly, the non-interacting bands (9.4) 
have the form: 



Ei{k) 
E2{k) 

E3{k) 

Ei{k) 



-E^{k) = -mv^ + E{k), 
+E-{k) = mv^ - E{k), 
+E+{k) ^ mv^ + E{k), 
-E+{k) = -mv^ - E{k), 



(9.5a) 
(9.5b) 
(9.5c) 
(9.5d) 



where E{k) = ^y {rnv'^Y + {vk)'^, and m — t±/{2v'^). 
Hence, Ei{k) and E^i^k) [or E2{k) and E^i^k)] describe 
a massive relativistic dispersion with rest energy given 
by mv^. Again, the gapless nature of the full spectrum 
of this problem is due to an accidental degeneracy of 
the simplest tight binding parametrization. Additional 
hopping terms (Castro Neto et ai, 2009a) in the Hamil- 
tonian or many body interactions can easily lift this de- 
generacy. This implies that the Bernal bilaycr problem is 
unstable from the electronic point of view. In contrast, 
the twisted bilaycr (Lopes dos Santos et ai, 2007) is sta- 
ble because it does not rely on this particular accidental 
degeneracy. Just like in the case of monolayer graphene, 
the introduction of the instantaneous Coulomb interac- 
tion does not preserve this Lorentz invariance. 

At very low energy, below A^^, w 1.5 meV, additional 
trigonal warping effects take place due to the influence 



of next-nearest neighbor hopping matrix elements [which 
we are neglecting in (9.2)]. Trigonal warping introduces 
an asymmetry in the conductivity under electron or hole 
doping (Li et ai, 2009b), and leads to a remarkable Lif- 
shitz transition at low densities, whereby the lowest en- 
ergy bands split into 4 Dirac cones (Cscrti et ai, 2007; 
McCann and Fal'ko. 2006). These effects, however, hap- 
pen at very low densities (around 1 electron per flake 
for typical samples), and hence are experimen- 

tally very challenging. A detailed description of the 
spectral properties of graphene bilayers can be found in 
Castro Neto et ai, 2009a, and Nilsson et ai, 2008. 

When Am < v\p\ <C t±, we recover the so-called clas- 
sical limit of the "relativistic" problem. This means that 
the presence of the uppermost band is not too relevant, 
and the energy disperses quadratically with momentum 
(the opposite limit of wlpj S> t_L corresponds to the 
"ultra-relativistic" regime, where the bandstructurc is es- 
sentially linear in momentum, like in the monolayer). In 
this case the Hamiltonian (9.2) near the K points can 
be projected onto an effective two-band model, written 
in terms of the two valleys and a mixed sublattice-layer 
basis (McCann and Fal'ko, 2006): 

^'p,ff = («K+p,(T, ^K+p,(T, &-K+p,(T, fl-K+p.o-) • (9.6) 

In such a basis, the effective kinetic Hamiltonian is 



pa a— ± 



■7) 



where p± = ipy, o± = (ui ± ?(T2)/2 operating in the 
sublattice basis, and r operates in the valley space. The 
resulting energy spectrum is parabolic. 



E{p) - ± 



2m 



(9.8) 



with m, = t±/{2v'^) w 0.054 TOg as the effective mass of 
the electron. From now on we will omit the valley indexes 
and assume the two component basis ^p.a (ap.o-i ^p.o-) 
with a total degeneracy = 4 in valley and spin. 

The electronic Green's function in this two band 
model, G(o)(k,T) = -(r[«'k(r)«'j^(0)]), is given by 
G^''^(k, iw) = {iuj — Hb)^^ or, equivalently, by 



G("'(k, 



-E- 



1 + S(Tk 



2 f-^^^u: - s\E{]^)\ 



in the chiral representation, where 



^2 

E iT^'^-- 



(9.9) 



(9.10) 



a=± 



Although the fermions are chiral, in bilayers the wave- 
functions of the quasiparticles acquire a 27r phase when 
winding around the K points, rather than a 7r-phase, as 
for Dirac fermions. This property is an admixture of the 
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FIG. 43 (Color online) The polarization n'-'^\q,uj) of bi- 
layer graphene, obtained within the two-band approxima- 
tion, for finite chemical potential, and zero temperature. All 
panels are normalized to the DOS at the Fermi energy, /i. 
Panel (a) shows a density plot of the imaginary part and, 
in (b), we have cuts of the same at constant frequency, for 
(jj/fj, — 0.5, 1.0, 1.5, 2.0, 2.5, 3.0. Panel (c) corresponds to the 
static limit 11*^^ (q, 0) in eq. (9.12), and includes the intra-band 
contribution (dashed), the inter-band contribution (dotted), 
and the full polarization (solid). In (d) we represent the real 
and imaginary parts of the polarization in the undoped case 
(9.14) as a function oi v — 2muj/q'^. 



behavior of Dirac particles, which are chiral, with con- 
ventional electrons, which disperse quadratically. The 
combination of chirality and a trivial Berry phase has 
a clear experimental signature in the suppression of the 
zero-level plateau in the quantum Hall effect of the bi- 
layer, whose plateaus are quantized by integer numbers 
(McCann and Fal'ko, 2006; Novoselov et ai, 2006). 



A. Charge polarization 

Within the two band model, the one loop polarization 
function has the generic form given in Eq. (2.12) for the 
single layer. The adaptations for the present case consist 
in considering the bilayer spectrum, and a new overlap 
factor, which, for the bilayer, reads 



'^.s,.s',p,q 



1, 



s'cos(26'p.p+q)]. 



(9.11) 



In this expression 0p,p+q is, again, the angle between the 
vectors p and p + q. Below we shall focus our discussion 
in terms of the effective two-band Hamiltonian (9.7), and 
dispersion (9.8). 



The polarization function 11*^"'^) (g, oj) at finite density 
was obtained by Hwang and Das Sarma, 2008a in the 
T = static limit. The full dynamical case was cal- 
culated by Sensarma et ai, 2010 at T = 0, and by 
Lv and Wan, 2010 at finite temperature. The finite den- 
sity result can be obtained in closed analytical form for 
T = 0; but, in order to avoid reproducing here those 
lengthy expressions, we simply present H'^-'((7,aj) graph- 
ically in Figs. 43(a,b). The explicit form of the static 
limit reads (Hwang and Das Sarma, 2008a; Lv and Wan, 
2010) 



H(i)((z,0) 



at zero temperature, with 



= .9(l^)-/(l^)%-2M (9-12) 



/(^) = ^v/^^ + ln(^^^^) (9.13a) 



2x 



9 



(x) = - V4 + a;4 - In 



(9.13b) 



The DOS at the Fermi energy, p(^[i) = Nni/{2TT), is con- 
stant and density independent, by virtue of the parabolic 
nature of the low energy approximations (9.7) and (9.8) 
[note, however, that the consideration of the full 4-band 
spectrum leads to a DOS which is linear in energy; in 
this sense, the correction to the DOS that arises from 
considering the 4 versus the 2 band model is not neg- 
ligible (Ando, 2007)]. In this sense the bilayer is simi- 
lar to the conventional 2DEG. However, just as in the 
monolayer, the existence of two symmetric bands adds 
an inter-band channel, leading to a rather different quasi- 
particle spectrum, in comparison with the 2DEG. This 
can be seen by directly comparing Figs. 4(b) and 43(a). 
The behavior of H'^^-'(g, 0) is shown in Fig. 43(c), together 
with its decomposition into intra- and inter-band contri- 
butions, which are respectively associated with the choice 
ss' = 1, or ss' = —1 in eq. (9.11). As intuitively ex- 
pected, the inter-band contribution dominates at large 
momenta/small densities, whereas the intra-band tran- 
sitions dominate the low momenta/large density regime. 
Unlike the monolayer, or the 2DEG, the polarization is 
constant for both q <^ kp and q ^ kp- The former limit 
makes the bilayer similar to the conventional 2DEG and 
monolayer graphene, while the latter is neither akin to 
the 2DEG (for which the polarization decreases rapidly 
with q/kp [Fig. 4(e)]), nor to the monolayer (for which 
it increases linearly [Fig. 3(e)]). Moreover, at precisely 
q — 2kp, H'^)(g,0) is sharply cusped, which contrasts 
with the behavior of a monolayer, whose derivative is 
continuous. According to the standard theories of linear 
response, this feature at 2kp has important implications 
for the behavior of the induced charge, the associated de- 
cay of the Friedel oscillations around charged impurities, 
the effective RKKY interaction among magnetic impuri- 
ties, Kohn's anomaly in the phonon dispersion, etc. For 
example, one expects qualitative differences between the 
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resistivity arising from Coulomb scattering in mono- and 
bilayer graphene: it should be stronger in the bilayer, 
and have a more pronounced temperature dependence 
(Hwang and Das Sarma, 2008a; Lv and Wan, 2010). 

At long wavelengths, the RPA screened potential reads 
V'^P^iq) = F(g)/[1-F(q)n(l)(q)] « 2^eV[eo(<Z + ^tf)] , 
with a Thomas- Fermi momentum qtf = Nme^ /tQ. No- 
tice that qxF is the same for the bilayer as in the 2DEG, 
i.e. it is constant (no density dependence), and also tem- 
perature independent (Lv and Wan, 2010). The temper- 
ature independence of qxp at long wavelengths is an- 
other trait that distinguishes this system from both the 
monolayer and the 2DEG. In real space the statically 
screened potential decays asymptotically as V{r) cx 1/r^ 
(Hwang and Das Sarma, 2008a). 

At half-filling (undoped situation) and zero tempera- 
ture, the form of the polarization bubble simplifies fur- 
ther, and can be cast as 



Nm 
"2^ 



v 



1 



2v 



l + 2v 
l-2v 



+ ln 



(9.14) 



(Barlas and Yang, 2009; Nandkishore and Levitov, 
2010a; Nilsson et ai, 2006), where v = 2ma;/g^ is the 
only scaling parameter. This function is plotted in 
Fig. 43(d). It follows at once that the static limit 
(w — >■ 0) is simply 



, TV In 4 

n(i)((?,o) = —m, 



(9.15) 



consistent with the above discussion when kp = 0. 
Despite the absence of a Fermi surface at half-filling, 
the Coulomb interaction among the quasiparticles is 
screened due to the finite density of states at the K 
points. However, an important difference here is that 
H(^^((jf,0) is constant for all momenta, unlike traditional 
2D systems, and stems from the presence of the inter- 
band channel. Hence, the Thomas-Fermi wavcvcctor 
is exactly qxp — iVm ln(4)e^/eo for all wavelengths, 
and Friedel oscillations are suppressed at half-filling 
(Hwang and Das Sarma, 2008a). The additional numer- 
ical factor ln(4) means a slight increase in the screen- 
ing strength of undoped bilayer, with respect to the 
doped situation. One way to interpret this ln(4) en- 
hancement is the following: the factor Nme^/eQ, being 
exactly the same as in a simple 2DEG, is attributable 
to the finite DOS, while the extra ln(4) arises from the 
virtual intcr-band transitions. In real space, the stat- 
ically screened potential of undoped bilayer will decay 
as which contrasts with the corresponding behav- 

ior in the monolayer, where the decay is 1/r (as we saw 
before this is due to the fact that, in the RPA, the ef- 
fect of interactions in the monolayer is to simply renor- 
malizc the background dielectric constant, keeping the 
Coulomb form of the potential). Inspection of Fig. 43(d) 



reveals that the real part of the RPA dielectric func- 
tion ei^p^(q, w) = eo[l — ^(q)n^^^ (q, w)] will be always 
nonzero. This means that, although the lack of a Fermi 
surface does not prevent screening in bilayers {qxp 7^ 0), 
the formation of zero temperature infrared plasmons is 
suppressed at half-filling. 

The screened Coulomb interaction between the lay- 
ers is V{q) = 27re^e~''''/[eo(g -I- q'tfC"^'')], where d — 
3.35 A is the interlayer distance. At long wavelengths, 
q <^ t_L/v < 1/d 0.3 A , d can be effectively replaced 
by zero in first approximation, and the screened interac- 
tion among electrons belonging to the same or different 
planes can be treated on the same footing. 

At this point we should pause to point out that the 
behaviors discussed so far at large q have to be inter- 
preted within the restrictions regarding the validity of 
the two-band approximation. For example, the fact that 
in Fig. 43(c) we see the polarization becoming constant 
at g 3> is an artifact of the two-band approximation. 
In reality, we should bear in mind that the full dispersion 
is hyperbolic, and hence becomes linear at high densities. 
We then expect to recover the linear-in-g dependence of 
n(^)(g,0) seen in Fig. 3(e) for the monolayer. 

For this reason, proper caution is needed when consid- 
ering the extrapolation of these results to highly doped 
bilayers, where the consideration of the four-band hy- 
perbolic dispersion (9.4) is inevitably required. In terms 
of electronic densities, this corresponds to values above 
^ 10^^ cm~^, for which the two-band model is no longer 
warranted. The full dynamical response using the spec- 
trum in eq. (9.4) has been recently derived in closed an- 
alytical form by Borghi et ai, 2009b. Notwithstanding 
the lengthy and cumbersome nature of these analytical 
results, they afford a more accurate perspective on the 
screening response of doped bilayer graphene, its col- 
lective modes, and the crossover between the regimes 
of a massive-chiral system at low densities, to a sys- 
tem of weekly coupled monolayers at higher densities. 
Borghi et ai, 2009b's approach is ultimately limited by 
systems of such high densities that /i « i, in which case 
the full tight-binding dispersion (9.3) is needed, but is 
beyond closed analytical approaches. 



B. Quasiparticles 

In the two band model, the structure of perturbation 
theory for Coulomb interactions is set only by self-energy 
renormalizations in the effective mass of the electrons, m, 
and in the quasiparticle residue, Z . 

From the Hamiltonian (9.7), the renormalized Green's 
function is 



G{k,Lo) = 



1 



Ea=± A:2/(2m)a„-I](k,a.) 



(9.16) 



S(k, w) is the quasiparticle self-energy correction, which 
is described in the (ak.o-, ^k.cr) basis by a matrix in the 
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2 4 6 8 10 

n (in units of lO"'^ cm"^) 



12 



for the two-band model (Borghi et al.. 2009a), and also 
for the four-band model (Kusminskiy et al., 2009), have 
found mass rcnormalization in the bilayer. The mass 
decreases (rn* /m < 1), and the renormalization grows 
stronger as the screening is suppressed. In Fig. 44 we 
show this rcnormalization within the two-band model, 
where the parameter A interpolates between the Thomas- 
Fermi screened potential (A = 1) and the unscreened 
Coulomb potential (A w 0). As a consequence of the re- 
duced mass, the charge compressibility is also expected to 
decrease (Borghi et al, 2010; Kusminskiy et al, 2008). 
More recent calculations that account for the 
screening have found quite differ- 
When the dynamical RPA polar- 
Eq (9.14), is taken into account, 



FIG. 44 (Color online) Mass renormalization for a = 0.5 in 
the bilayer, calculated with a static Thomas- Fermi screened 



full dynamical 
ent results, 
ization bubble, 
the self-energy 
gence, KeYS^\k.uj 
(1) 



exhibits a strong In leading divcr- 
= 2fc2/(iVm7r2)ln^(A//s), and 



of the electronic density. Blue circles; A = 1; red squares: A 
0.01; green triangles: A = lO""*. (From Borghi et al, 2009a). 
The inset shows logj^(j(7Ti*/m) as a function of logjQ(n) for two 
of the A values; the mass saturates at a finite value for n — >■ 0. 



form 



Coulomb interaction V{q)^= e /[eo{q + Mtf)], as a function ReSg ' [k,uj) = -Alo / {NTT'^)\n [K/ {^/rnuj)], at smaU en- 
ergies and momenta (Barlas and Yang, 2009). The ultra- 
violet momentum scale A ^ q^F is related to the effective 
"Bohr radius", ao = eo/ime?), and we set A = l/flQ. At 
leading (In^) order, the two terms in the self-energy com- 
pensate each other exactly in Eq. (9.19) and the mass 
does not renormalize, m* /m — >■ 1 at fc — > 0, while the 
quasiparticle spectral weight vanishes as Z ~ ln~^(A//c). 
The RG analysis of the dynamically screened interaction 
at large N was carried out by Nandkishore and Levitov, 
2010c, where subleading (single log) contributions were 
collected. These were found to cause a (weak) increase 
of the effective mass m* /m sa 1 -I- [0.56/(A^27rln4)] In A, 
and consequently an increase of the compressibility. 

Once again, the validity of a two band model rests 
on the assumption that all relevant energy scales are 
small compared to t± ~ 0.4eV. However the Coulomb 
energy A^; on the scale of ao = eo/(me^) is substantial 



So S+ 
S_ So 



(9.17) 



or, equivalently, E = Sqcto + 5]+<t+ + S_(t_, where 
a± = (cr^ ± iay)/2. By symmetry, = El. In a more 
conventional form. 



G'(k,w) = 



^Ea=±[fc^/(2m) + EJa<,' 



(9.18) 



where Z ^ = 1— SEq/Sw corresponds to the quasiparticle 
residue, and 



m* _ 1 - dT.n/duj 
~^ ^ 1 + 2mdY.+/dk\ 



(9.19) 



is the mass renormalization. 

We saw in the previous section that, unlike the mono- 
layer. Coulomb interactions in the bilayer are screened. 
The self-energy is given in terms of the bare Green's func- 
tion and the RPA effective interaction by 

cPkde 



for not too strong dielectric screening, A^; = e /(eooo) ~ 
1.47/e§ eV (Nandkishore and Levitov, 2010a,c). Hence, 
Coulomb interactions can promote electronic transitions 
among the four bands, while the two-band model is only 
justified in the limit Ke < t^. To what extent the two 
band model provides a valid description of the quasipar- 
ticlcs in the presence of Coulomb interactions is a matter 
of ongoing discussion. 



\ ■ f ^rRPAn \/^(0)/^ , , \ C. Many-body instabilities 

E^'^ (q, w) I y -^^-y3-y"^^(k, e)G^''^ (k q, e w), 



(2^) 

(9.20) 

where V^^^{q,uj) = F(q)/e/jpA(q, i^) is dressed by the 
RPA dielectric function. Even if the ratio between the 
Coulomb and kinetic energies diverges in the low density 
limit (as in a 2DEG) the validity of RPA can be, in prin- 
ciple, justified in the large N limit. If only static screen- 
ing is taken into account (Hartree-Fock- Thomas-Fermi 
theory), the self-energy is frequency independent and, to 
leading order, the quasiparticle residue Z does not renor- 
malize. Calculations based on the static screening picture 



The finite DOS in the bilayer enhances the possibility 
of many body instabilities in comparison with the single 
layer case. For instance, the spin polarization tensor in 
the bilayer is defined in leading order by Eq. (5.4), with 

the matrix element ^^(k) = 1 + sJ2a=±(^a/^^)'^a- 
matrix form. 



n(i)- 



H 



(1) 



n^''+" H 



b 

(1) + - 



(9.21) 
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which leads to one ferromagnetic, and one antiferro- 
magnetic eigenstate, Uf/af = ^aa ± IIIq;,!, by sym- 
metry under exchange of the a and b labels. In bilay- 
ers the AF state has a leading logarithmic divergence 
with the cut-off, A, at zero frequency and magnetic field, 
(Nilsson et al, 2006) 

nW(q,0) = ^ln(^), (9.22) 

suggesting (within RPA) a tendency towards an AF in- 
stability for any value of the Hubbard interaction U . In 
addition, at finite U , a first order ferromagnetic tran- 
sition can be driven by the Stoner criterion, leading to 
a ferromagneto-electric state where the layers have dif- 
ferent magnetization and polarized charge (Castro et al., 
2008b). 

Other possibilities include the emergence of CDW 
instabilities induced by the short range part of the 
Coulomb interaction, (Dahal et al., 2010) or else, an 
excitonic instability at strong local electronic repul- 
sion (Dillenschneider and Han, 2008). With long 
range Coulomb interactions, the inverse electronic com- 
pressibility becomes negative at small densities 
(Kusminskiy et al., 2008), indicating a tendency to 
Wigner crystallization (Dahal et al., 2006), which is com- 
pensated by the positive compressibility of the lattice. 

Bilaycrs share similar features with one-dimensional 
(ID) electron systems, such as the point like Fermi sur- 
faces and the parabolic spectrum. In particular, in bi- 
ased bilayers, the ID interface between biased regions 
confines chiral modes that propagate as in a strongly in- 
teracting Luttinger liquid (KiUi et al., 2010). This af- 
fords the possibility of studying such interacting models 
experimentally in appropriately prepared samples of bi- 
layer graphcne. 

For short-ranged interactions in 2D, the structure of 
the diagrams in bilayers and in ID electron liquids is quite 
similar, although the diagrams compensate each other in 
a rather different way. The dimensionless coupling which 
determines the strength of the interactions is Ua^m, 
where U is the strength of the local interactions, and a is 
the lattice constant. Pcrturbativc rcnormalization group 
calculations in the bilayer have identified distinct leading 
instabilities of the electron gas. For different choices of 
possible interactions, two different low-temperature bro- 
ken symmetry phases have been found: in one case, a fer- 
roelectric gapped phase (Zhang et al., 2010) induced by 
the coupling between the different layers; in the other, 
a nematic phase (Vafek, 2010; Vafek and Yang, 2010), 
where each Fermi point splits into two Dirac points. 

The possibility of an excitonic instability has been also 
predicted by Nandkishorc and Levitov, 2010a, who found 
that the dynamically screened Coulomb interaction gives 
rise to a ferroelectric state that polarizes the two lay- 
ers. In the ferroelectric state, the kinetic energy inflicts 
an energy cost <5ii^Kinetic oc ln(A£;/A), where A is the 
energy gap. Finite separation between the layers gener- 
ates an additional electrostatic energy cost to polarize the 



charge between the layers, which dominates the kinetic 
energy at the Hartree level, (J-Enartrco c>c A^ln^(A£;/A) 
(McCann et al., 2007). The excitonic instability is in- 
duced by the exchange term, which is parametrically 
larger than the Hartree term by the factor (ao/d), where 
d is the interlaycr distance (Nandkishorc and Levitov, 
2010a). The existence of a ferroelectric state has nev- 
ertheless been disputed by independent RG calcula- 
tions, that also accounted for the dynamically screened 
Coulomb interactions, and infrared trigonal warping ef- 
fects (Lcmonik et al., 2010). The spontaneous sym- 
metry breaking found in this work leads to a Lifshitz 
transition consistent with the nematic state found by 
Vafek and Yang, 2010, rather than the opening of a gap. 

In the Quantum Hall (QH) state, two terminal mea- 
surements of the conductivity in clean suspended sam- 
ples have found an insulating state at the v = fill- 
ing factor (Feldman et al., 2009), rather than the metal- 
lic QH state previously found in supported samples 
(Novoselov et al, 2006). Further theoretical works pre- 
dicted the possibility of a zero field excitonic QH state, 
which spontaneously breaks time reversal symmetry, and 
can evolve into a ferromagnetic QH state at finite mag- 
netic field (Nandkishorc and Levitov, 2010b). In biased 
bilayers, a chiral anomaly has been predicted in the 
QHE, splitting the degeneracy of valley quantum num- 
bers (Nakamura et al., 2009). Another predicted effect 
resulting from interactions in the QH state is the for- 
mation of charge 2e skyrmions at even filling factors 
(Abanin et al, 2009). 



X. CONCLUSIONS 

As wc have seen, the understanding of the many-body 
problem in graphcne has evolved quite rapidly in only a 
few years. The case of monolayer graphene in the weak 
coupling regime (which means, graphene embedded in 
an environment with large dielectric constant) is quite 
clear, namely, although Lorentz invariance is explicitly 
broken because of the Coulomb interactions, the effec- 
tive low energy theory is still Lorentz invariant with well 
defined quasiparticles. Nevertheless, these quasiparticles 
have a renormalizcd speed of light that grows logarithmi- 
cally in the infrared, while their spectral weight decreases 
slowly in the same limit. This situation can be contrasted 
with the conventional Fermi liquid picture where all the 
physical constants (the so-called Landau parameters) and 
spectral weight are finite in the infrared (that is, at the 
Fermi surface). Hence, these logarithmic renormaliza- 
tions are weak enough, even in the presence of strong 
Coulomb interactions, and a Dirac liquid picture is pre- 
served. 

In the strong coupling regime (that is, graphene in 
vacuum), many-body instabilities are possible albeit de- 
pending on a delicate balance of energy scales. This oc- 
curs because the renormalizations of quasiparticle prop- 
erties also depend on details of the cut-off procedure in 
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the ultraviolet (as it is shown by the f-sum rule). While 
mean-field theories have predicted instabilities towards 
phases with broken chiral symmetry and superconduct- 
ing quasi-long range order (because of the 2D nature 
of the material), and earlier Monte Carlo studies on a 
hyper-cubic lattice suggest the presence of instabilities 
(Drut and Lahde, 2009a,b,c), simulations of interacting 
electrons on the honeycomb lattice have still to be per- 
formed in order to address these issues, since the strong 
coupling regime cannot be reached by perturbative meth- 
ods. This remains, at this point in time, as an important 
open problem in many-body graphene physics. 

The Coulomb impurity problem in graphene shares 
many of the issues of the many-body problem but can 
be studied in much more detail because the 2D hydro- 
gen problem in graphene was solved exactly. In the weak 
coupling regime (the so-called under-critical regime), the 
Coulomb interaction between a localized charge and the 
electrons leads to only mild changes in the physical prop- 
erties due to the explicitly broken particle-hole symme- 
try. In the strong coupling (or super-critical) regime, the 
situation is rather different because of the phenomenon of 
fall to the center, that is, the electron states become un- 
stable, with the generation of resonances near the Dirac 
point. Just like the many-body problem, the critical lo- 
cal charge depends on the dielectric environment, and in 
vacuum this amazing effect should be observed by local 
probes even for a single proton sitting on the graphene 
surface. So far, there is no experimental evidence of such 
effect, given that it is difficult to study adatoms in sus- 
pended samples with local probes, such as scanning tun- 
neling microscopes. In supported samples, because of 
dielectric screening that brings the system to weak cou- 
pling, and of the disorder in the substrate, the study of 
this problem can be much more elusive. 

In analogy to the 2DEG problem, the effect of disor- 
der is rather strong in graphene which again is the ef- 
fect of dimensionality. The low dimensionality implies 
strong quantum fluctuations that can easily couple to 
spatial variations of random scalar (chemical potential) 
and gauge (hopping) fields. Strong localization is the 
ultimate fate of any disordered two dimensional system 
but because the localization length grows very slowly in 
the infrared limit, the finite size of the samples, or the 
finite temperature of the system, ends up cutting off the 
tendency towards Anderson localization and, in practical 
terms, graphene behaves in a metallic way. 

The problem of magnetism of adatoms in graphene is 
rather different from the one found in metallic hosts. Due 
to the strong energy dependence of the density of states 
(that vanishes at the Dirac point), the Anderson impu- 
rity problem has features that are unique. Firstly, in 
analogy with the strong coupling regime in the many- 
body and Coulomb impurity problems, the results are 
sensitive to the ultraviolet regularization. In fact, this is 
a generic feature of the Dirac spectrum, namely, strong 
coupling leads to spectral weight transfer from high en- 
ergies to low energies, that is, to the Dirac point (as 



discussed in the context of the f-sum rule). Moreover, 
the damping by Dirac electrons leads to an anomalously 
large (and strongly energy dependent) broadening of the 
adatom energy level. This leads to an unusual situa- 
tion as compared to the Anderson impurity problem in a 
metal, namely, that even when the chemical potential is 
above (below) the energy of the doubly (singly) occupied 
state, a magnetic moment can emerge. Hence, adatoms 
that may not be magnetic in a metal (hydrogen or fluo- 
rine, for instance), might become magnetic in graphene. 

On the other hand, the Kondo effect that usually sup- 
presses the appearance of magnetic moments in metals 
because of magnetic "screening" (the ultimate conse- 
quence of the so-called "Kondo cloud" ) , is strongly sup- 
pressed in graphene. This suppression has its roots in 
the low density of states and the sublattice structure. In 
fact, there is a strong dependence of the hybridization 
with the position on the lattice (whether it breaks or not 
the sublattice symmetry) . Furthermore, the Kondo effect 
is very dependent on the chemical potential (that can be 
easily tuned in graphene by gating) . This state of affairs 
reinforces the conclusion that magnetic states of adatoms 
could be more the norm than the exception in graphene, 
in a big contrast with the situation in ordinary metals. 
Experimentally, there are very few studies of the mag- 
netism of adatoms in graphene. The main problem here 
is that most of the experiments done so far are in elec- 
tronic transport. Just like the Kondo problem in metals 
and semiconductors, the observation of magnetic effects 
in transport is rather subtle, and requires careful analy- 
sis. At this point in time, this is a rather open field in 
graphene physics. 

A superconducting state in graphene would have dra- 
matic consequences given its low dimensionality and un- 
usual electronic spectrum. While true long range order 
would not be possible because of its 2D nature, quasi- 
long range order would have unusual consequences. For 
one, because of the sublattice structure, there is room for 
exotic pairing states with even more exotic vortex excita- 
tions. The phase space for pairing is rather large due to 
the spin, sublattice, and valley degeneracies. However, 
the low density of states plays a deleterious role here. 
One way out of this conundrum would be the enhance- 
ment of the density of states by either gating or doping 
with adatoms. These two techniques have their own lim- 
itations. Gating is limited by the distance from the gate 
to the graphene sample, and by the dielectric breakdown 
of the spacer that separates the two. Doping inevitably 
introduces disorder, or can modify the electronic struc- 
ture of the TT band too much leading to extrinsic effects. 
There are, however, serious hopes that come from the fact 
that intercalated graphite can be made to superconduct. 
An obvious idea would be intercalation of Ca or Yb in 
the graphene bilayer. So far, intercalation experiments 
in bilayers have not been performed, and very little is 
known about how to intercalate atoms or molecules in 
such systems. Again, this is very much an open field of 
research. 
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In addition to the dielectric environment, which has a 
strong influence on many-body effects in bulk graphene, 
finite size effects are also of great importance. It has been 
understood very early on that zig-zag edges are strongly 
interacting because of the high density of states they cre- 
ate at the Dirac point. Systems with high density of 
states are prone to many-body states due to Stoncr-like 
instabilities. However, the many-body physics of finite 
2D systems is even more sensitive to disorder (either in 
the bulk or in the edge) because of the strong bound- 
ary condition dependence. In graphene this problem is 
magnified because the electronic wavefunctions associ- 
ated with impurity states do not decay exponentially, as 
they would in a semiconductor with a finite gap, or would 
be extended, as in a normal metal, but they are quasi- 
localized (that is, decay like power law). This implies 
that evanescent waves play an important role in deter- 
mining the physical properties. Experiments in meso- 
scopic graphene samples show very clearly these effects 
through strong oscillations of the electronic conductance 
and the presence of Coulomb blockade peaks. From theo- 
retical perspective, such problems arc probably the hard- 
est to solve because they involve the direct interplay 
between Anderson localization and interactions. Thus 
deeper understanding of mesoscopic graphene systems is 
still necessary, and this topic would merit a review of its 
own. 

Magnetic fields also lead to spatial localization due to 
the presence of Landau levels with a length scale given 
by the cyclotron length. Hence, this problem shares 
many of the difficulties of the previous problems with 
the added complication that the 2D nature of graphene 
brings a huge degeneracy into play. Once again, the 
detailed balance between kinetic and Coulomb energies, 
and the details in the ultraviolet, determine the fate of 
the many-body ground state. The fractional quantum 
Hall effect was only observed recently in suspended two- 
probe experiments (Du et al, 2009), and very little is 
known about the sequence of FQHE fractions and their 
nature. It is believed that magnetic fields can generate a 
plethora of new many-body states, with symmetries that 
arc rather different from the ones found in the 2DEG. 
But, compared to the 2DEG problem, this field is still in 
its infancy. 

While we have demonstrated the complexity of the 
many-body problem in monolayer graphene, we have not 
even touched beyond the surface of the many-body prob- 
lem in bilayer graphene. There is no doubt, at least from 
the theoretical perspective, that the many-body problem 
in the bilayer is much richer than in the monolayer. For 
one, the bilayer has a finite density of states at neutral- 
ity, making it similar to the 2DEG problem. However, 
unlike the 2DEG, the graphene bilayer is a Lorentz in- 
variant system with a finite "rest mass" (that is, it has a 
hyperbolic dispersion relation) albeit with an accidental 
degeneracy that makes it a semi-metal (two of the four 
bands touch at the Dirac point). This accidental degen- 
eracy can be lifted very easily by hopping or interactions. 



leading to a huge number of possible many-body states 
with different quantum numbers. Given this richness one 
can venture saying that bilayer graphene is the ultimate 
target of many-body theorists in this field. However, it is 
technically a major challenge given the high dimension- 
ality of the problem, with its 2'*-dimensional spinorial 
structure (spin, valley, sublatticc, and plane). Moreover, 
from the experimental perspective many details and con- 
ditions are still quite uncontrolled, which has led to a 
few contradictory results, and has so far yielded more 
questions than answers. In fact, both theoretically and 
experimentally, the graphene bilayer remains very much 
an open problem. If we now extrapolate from the mono- 
layer to the bilayer, we see that there are problems that 
have not even been addressed theoretically and exper- 
imentally, like the Anderson impurity problem, or the 
Kondo effect in bilayers, the problem of magnetism, and 
superconductivity, just to mention some. These are top- 
ics for the future, for future generations of physicists to 
address and marvel. 
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